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ABSTRACT 

This  paper  describes  a user-friendly  computer  model,  TMPSUB2.  This  model  calculates  the  temperature 
field  throughout  a solid  (but  porous)  substrate  which  can  undergo  exothermic  pyrolysis,  when  it  is 
exposed  to  an  arbitrary  (localized)  heat  flux  which  can  move  uniformly,  such  as  that  of  a smoldering 
cigarette.  TMPSUB2  has  successfully  simulated  the  thermal  runaway  signifying  smoldering  ignition  of 
the  substrate  when  it  is  exposed  to  a set  of  external  heating  fluxes.  The  processes  taken  into  consideration 
are  three-dimensional  heat  conduction  in  the  substrate,  pyrolysis  of  the  latter,  and  the  diffusion  of  air  into 
it  at  the  top  surface.  TMPSUB2  also  takes  into  account  the  fact  that  the  substrate  consists  of  two  layers 
(a  foam  pad  covered  by  a fabric),  with  the  possibility  of  an  air  gap  between  them,  up  to  three  pyrolytic 
reaction  steps,  and  with  temperature-dependent  thermophysical  constants.  TMPSUB2  solves  the  equations 
describing  the  physics  and  chemistry  of  the  heating  and  ignition  process  numerically;  the  results  have 
been  compared  with  a set  of  ignition  experiments,  and  have  been  found  to  be  semi-quantitatively  correct, 
both  for  the  ignition  temperature  and  for  the  time  to  ignition.  Analysis  of  the  experiments  indicates  that 
the  substrate,  which  consists  of  a thin  fabric  layer  over  a thick  foam  padding,  behaves  as  a thermally  thin 
layer.  Use  of  TMPSUB2  shows  that  smoldering  ignition  would  result  from  a stationary  hot  spot  of 
intensity  and  dimension  simulating  a quietly  smoldering  cigarette.  A users’  guide  is  included. 


Key  Words:  Computer  models;  furniture  fires;  ignition;  mathematical  modeling;  pyrolysis;  smoldering 


1.  INTRODUCTION 

Over  1600  people  die  each  year  (see  ref.[l])  from  fifes  started  by  smoking  materials  (cigarettes,  matches, 
etc.),  and  many  times  that  number  are  hurt  and/or  disfigured;  most  of  those  fires  are  established  in 
upholstered  furniture.  As  a result.  Congress  passed  the  Cigarette  Safety  Act  of  1984;  its  goal  was  to 
determine  whether  safer  (from  the  point  of  view  of  fires!)  cigarettes  were  technically  feasible.  One  of 
the  tasks  under  that  Act  was  to  model  mathematically  the  behavior  of  upholstered  furniture  and  of  the 
cigarette,  when  a smoldering  cigarette  is  dropped  onto  it.  The  result  was  two  prototype  programs, 
TEMPSUB  (for  TEMPerature  of  a SUBstrate)  and  CIG25  (see  ref  [2]).  One  of  the  tasks  under  the  Fire 
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Safe  Cigarette  Act  of  1990  (P.L.  101-352)  is  to  "Conduct  laboratory  studies  on  and  computer  modeling 
of  ignition  physics  to  develop  valid,  "user-friendly"  predictive  capability."  In  order  to  do  this,  two  (new) 
models  were  to  be  developed:  one  for  the  ignition  of  the  substrate  (an  advance  over  TEMPSUB),  and 
another  for  the  quiescent  smoldering  of  a cigarette  which  is  resting  on  the  substrate  (an  improvement  on 
CIG25).  In  this  paper,  the  first  model  is  described;  the  cigarette  model  will  follow. 

Consider  what  is  involved  in  developing  this  computer  model:  we  must  simulate  the  behavior  of  a 
(typical)  substrate  when  it  is  subjected  to  a heat  flux.  In  order  to  do  that,  the  physics  and  chemistry  of 
the  process  must  first  be  understood.  It  is  then  expressed  as  a set  of  equations  describing  the  behavior; 
these  equations  must  then  be  solved,  with  the  appropriate  data  and  boundary  conditions.  The  solution 
method  is  numerical,  and  the  data  input  is  to  a computer. 

This  program  serves  to  calculate  the  temperature  of  the  upholstered  furniture  as  a function  of  time  and 
position,  when  it  is  exposed  to  a prescribed  heating  flux.  This  flux  can  be  highly  peaked  at  a point,  vary 
with  time,  and  move  at  a constant  (specified)  rate  over  the  top  surface  of  the  furniture,  assumed  to  be 
horizontal.  The  radiative  and  convective  heat  losses  from  the  surface  are  given  correctly.  If  and  when 
the  temperature  rate  of  rise  at  a given  location  suddenly  "accelerates"  to  a value  high  enough  that  the 
surface  glows  (that  is,  T > 5(X)°C  or  so),  we  can  say  that  smoldering  ignition  has  occurred.  The 
ambient  oxygen  level  can  be  set  at  whatever  value  one  wishes.  The  program  will  not  tell  us  whether 
flaming  ignition  takes  place.  It  also  does  not  treat  the  case  where  the  flux  is  applied  in  a crevice,  such 
as  is  formed  between  the  seat  cushion  and  the  seat  back.  The  program  does  not  take  oxygen  diffusion 
within  the  cushion  explicitly  into  account;  hence  in  certain  threshold  situations,  where  a small  change  in 
oxygen  concentration  determines  whether  ignition  does  or  does  not  take  place,  the  results  are  ambiguous 
and  not  to  be  trusted.  Note  that  it  is  often  difficult  to  obtain  the  needed  kinetic  and/or  thermophysical 
parameters  for  the  material;  or,  when  available,  to  know  how  accurate  they  are.  Therefore  this  caveat 
must  also  be  made:  even  if  the  program  were  perfect,  its  results  are  only  as  good  as  the  input  parameters 
which  are  supplied.  On  the  other  hand,  it  should  accurately  reproduce  (or  predict)  trends. 

Here  is  a brief  overview  of  the  dynamics:  For  the  purpose  of  this  work,  the  "upholstered  furniture"  will 
be  taken  to  be  the  flat,  horizontal  cushioned  seat  - that  is,  the  substrate  is  a cushion  consisting  of  fabric- 
covered  foam  padding.  Assume  that  there  is  a localized  source  of  heat  on  the  cushion.  The  temperature 
of  the  cushion  will  then  rise,  locally.  In  order  for  ignition  to  occur,  a number  of  conditions  must  be  met; 
the  principal  one  is  that  the  temperature  of  the  cushion  must  rise,  in  this  region,  to  a critical  temperature 
for  the  initiation  of  smoldering.  This  is  generally  referred  to  as  the  ignition  temperature.  It  will  be 
shown,  however,  that  this  temperature  is  not  unique,  even  for  a well-defined  material:  it  depends  also 
on  the  magnitude  of  the  incident  (igniting)  flux.  The  temperature  of  the  substrate  will  be  calculated, 
nevertheless;  for,  as  will  be  shown,  the  temperature  history  will  determine  whether  ignition  does  or  does 
not  take  place. 

In  order  to  calculate  the  temperature  as  a function  of  time  and  position,  it  is  necessary  to  know  the  heat 
flux  from  the  cigarette  (as  it  is  affected  by  its  contact  with  the  chair),  how  this  heat  flux  distribution 
moves  with  time,  and  how  the  seat  cushion  responds  to  that  "insult."  We  wish  to  find  whether  the 
cushion  will  in  fact  reach  its  ignition  temperature,  to  begin  with.  It  may  be  necessary  to  know  some 
other  things  as  well,  such  as  the  effect  of  the  presence  of  the  cigarette  on  the  oxygen  concentration  in  the 
cushion;  that  will  be  dealt  with  at  a later  time.  Thus  the  first  task  is  to  examine  the  thermal  response  of 
the  cushion  to  a specified  flux  distribution.  This  is  done  in  the  next  section.  In  section  3,  the  numerical 
solution  method  is  presented,  and  some  of  the  physics  is  revisited  in  this  (numerical)  context.  In  section 
4 we  discuss  a set  of  experiments  where  a mock-up  was  ignited,  give  the  measured  ignition  times,  and 
give  the  input  data  requited  to  make  computer  runs.  The  results  of  these  runs  are  given  and  then 
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discussed  in  Section  5,  validating  the  model  to  some  degree.  Section  6 summarizes  the  work.  For  the 
reader  interested  in  a more  extensive  discussion,  see  the  references  cited  here,  especially  the  bibliography 
in  [2].  A copy  of  the  program  on  disk  may  be  obtained  from  the  authors,  as  well  as  a listing  of  the 
program  and  the  program  itself. 


2.  IGNITION  DYNAMICS 

The  equations  describing  the  relevant  physical  processes  taking  place  in  the  substrate  will  be  written  in 
this  section.  They  will  mostly  be  written  in  general  form  to  start,  and  then  particularized  to  our  problem. 
The  assumptions  made  in  the  model  are  summarized  and  displayed  in  Table  1 . 

Table  1.  What  is,  and  is  not,  in  the  substrate  model 


IN 

• two  porous  layers  with  an  air  gap 

• 3-d  heat  conduction 

• variable  (prescribed)  thermophysical  properties 

• correct  (nonlinear)  boundary  conditions 

• endothermic  (non-oxidative)  pyrolysis 

• exothermic  (oxidative)  pyrolysis 

• char  oxidation 

• arbitrary  (prescribed)  moving  heat  source 

• variable  grid  on  a moving  coordinate  system 

• (approximate)  radiative  heat  transfer  within  the  material 

• Impinging  air  flow,  when  wanted 

NOT  IN 

• oxygen  diffusion 

• melting  and/or  regression  of  the  foam 


The  most  appropriate  place  to  begin  (perhaps)  is  with  the  heat  source;  that  will  be  described  in 
subsection  2a.  The  equation  describing  heat  diffusion  in  a solid  is  given  and  discussed  in  subsection  2b; 
it  is  solved  for  a very  simple,  special  case  in  section  2c,  for  later  use.  Pyrolysis  of  the  substrate  is 
discussed  in  section  2d,  and  the  diffusion  of  gases  (a  process  very  similar  to  the  diffusion  of  heat)  in 
section  2e.  The  principal  gas  of  interest  here  is  oxygen.  Although  we  do  not  explicitly  include  its 
movement  in  the  substrate,  those  expressions  are  necessary  for  estimating  the  rate  at  which  oxygen  can 
enter  the  reactive  zone  in  the  substrate,  from  outside.  That  boundary  condition  and  its  consequences  is 
discussed  in  subsection  2f. 

2a.  Surface  Heat  Source 

The  process  is  initiated  by  external  heating  of  the  surface.  Thus,  it  is  necessary  to  specify  the  heat  gains 
at  the  bounding  surfaces  of  the  cushion.  The  cushion  is  assumed  to  be  a rectangular  parallelepiped; 
therefore  a Cartesian  coordinate  system  is  employed.  For  the  top  surface,  we  have:  first,  the  localized 
heat  flux  from  the  glowing  tip  of  the  cigarette  or  other  heat  source;  the  flux  consists  of  two  parts,  one 
due  to  convection  (<t>^  and  one  due  to  radiation 
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(1) 


i>i„  = 4),  + 4), 

is  the  heating  flux  reaching  the  surface.  For  a cigarette,  these  fluxes  are  strong  functions  of  position, 
since  the  glowing  tip  is  only  a few  millimeters  in  extent.  With  this  formulation,  <f)^  can  be  expressed  in 
the  approximate  form 

4.,  - Q€^,»r4-(i-o)or/  C) 


where  = Q(r)  is  the  shape,  or  view,  factor  of  the  cigarette  as  seen  by  the  substrate  at  the  point  r.  It 
is  only  approximate  in  this  form,  because  the  cigarette  surface  temperature,  Tjg,  varies  strongly  with 
position. 

The  convective  flux  is  given  by 

4>,(^,y.O  = h[T^ig(x,yyt)  - T/x,y,t)]  (3) 

where  h is  the  heat  transfer  coefficient,  is  the  surface  temperature  of  ±e  substrate,  which  is  explicitly 
shown  to  be  a function  of  position  and  time,  and  where  Tj,jg  is  the  cigarette  surface  temperature,  also  a 
function  of  position  and  time.  It  is  important  to  point  out  here  that  in  this  version  of  the  substrate 
model,  the  source  flux  must  be  specified  by  the  user. 

It  is  equally  necessary  to  specify  the  heat  losses.  If  the  existence  of  the  cigarette  is  neglected  for  the 
moment,  so  that  the  view  factors  need  not  be  considered,  then 

= (4) 


where  the  first  term  is  the  convective  cooling  term  and  the  second  is  the  radiative  cooling  term.  In 
contrast  to  the  source  term,  these  losses  are  explicitly  included  in  the  program.  Here,  the  dependence 
of  Tg  on  X,  y,  and  t has  not  been  shown  explicitly,  e is  the  emissivity  of  the  surface  and  o is  the  Stefan- 
Boltzmann  constant  (note:  we  have  used  e where  sometimes  a,  the  absorptivity  of  the  surface,  should 
be  used.  When  there  is  thermal  equilibrium  between  the  radiation  field  and  the  hot  surface,  then  a = 
e.  We  have  assumed  the  latter,  for  simplicity).  Evidently,  the  loss  rate  from  the  cushion  surface 
increases  as  it  heats  up.  h is  determined  by  the  laws  of  fluid  flow,  as  well  as  the  thermal  properties  of 
air;  the  simplest  way  to  find  h,  however,  is  to  use  well-known  expressions,  usually  derived  from 
correlations.  This  will  be  further  discussed  below. 


Along  with  the  partial  differential  equation  (12),  we  must  have  boundary  conditions:  for  the  top  surface, 
the  net  flux  entering  the  surface  at  each  point  is  connected  to  the  net  flux  according  to 


dT{x,y,z,t) 


dz 


z = 0 , t>0 


(5) 


z = 0 


where 


4*net  4*in  “ 4*out  • 


(6) 


When  the  cigarette  lies  on  the  substrate,  the  mean  convective  heating  flux  over  the  heating  region  is  given 
by 


4>  = A- 

II 


(T  . -T) 

r CIg  s' 


(7) 
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If  hq  is  the  heat  transfer  coefficient  for  quiescent  air,  then  the  surface  heat  loss  for  the  areas  away  from 
the  cigarette  is 

d)  = /i  (7  - r ) (8) 

^out  s o' 

The  boundary  condition  is  greatly  simplified  if  we  can  say  that  Eq(8)  is  valid  over  the  entire  surface. 
That  is  readily  achieved  by  using  the  model  convective  heating  flux  <!>*, 

d)’  = A.  (f  . - T)  + /i  (7  - r ) (9) 

▼c  eig  s'  9''  * o' 

hjn  is  found  as  follows:  the  mean  initial  convective  flux  from  a cigarette  was  measured  to  be  37  kW/m^. 
Taking  the  mean  value  of  the  cigarette  surface  temperature  to  be  about  <T>cjg  = 450°C,  the  mean  heat 
transfer  coefficient  in  this  area  is  h^^  = (37000/430)  = 86  W/m^-K.  h^  is  found  from  Table  7-1 
of  ref  [2a]:  the  Nusselt  number  for  a heated,  horizontal,  upward-facing  surface  is 

Nu  = OMiRay* 

where  Ra  is  the  Rayleigh  number.  Since  h = kNu/L,  where  L is  the  characteristic  dimension  of  the 
surface,  we  find 


h 


0.54  k 


g^iT-TJ 


avL 


V4 


(11) 


where  jS  is  the  volumetric  coefficient  of  (gas)  expansion.  Taking  the  value  for  T^  = 500  K = 227  °C 
as  an  approximate  mean  temperature  of  the  surface,  and  L = 7 mm,  we  find 


\ = 9.7  W/m^-K 

For  the  other  surfaces,  the  boundary  conditions  (b.c.)  can  be  expressed  in  different  forms,  depending  on 
whether  the  slab  for  which  we  are  making  calculations  is  exposed  to  the  air,  or  is  embedded  within  the 
cushion.  If  the  slab  were  the  entire  cushion  - that  is,  the  other  five  surfaces  are  exposed  to  the  air,  then 
the  b.c.  would  be  that  the  heat  fluxes  KdTIdx,  xdT/dy,  KdT/dz,  are  given  by  terms  of  exactly  the  same 
form  as  (4),  except  that  h is  different  for  the  vertical  and  the  downward-facing  horizontal  surface.  On 
the  other  hand,  the  numerical  calculation  we  are  using  makes  using  the  entire  cushion  prohibitively  large. 
We  therefore  consider  a subsection  abstracted  from  the  whole  cushion,  and  the  other  five  faces  are 
surrounded  by  more  foam.  If  that  foam  were  a perfectly  insulating  material,  then  the  b.c.  at  those  faces 
is  the  adiabatic  condition;  that  is,  no  flux  crosses  the  surfaces:  xdT/dx  = xdTIdy  = KdT/dz  = 0.  This 
eliminates  any  heat  losses  from  the  slab,  and  makes  the  calculated  temperatures  somewhat  higher  than 
they  really  are.  A simpler  b.c.  is  to  assume  that  the  temperature  at  the  five  faces  is  constant  and 
unchanging,  remaining  at  the  original  (usually  ambient)  temperature.  This  produces  a larger  temperature 
gradient  than  actually  prevails,  and  errs  in  the  opposite  direction.  That  is  the  default  b.c.  used  in  the 
calculation,  but  which  of  the  two  b.c.’s  is  used  is  determined  by  an  input  parameter. 


2b.  Heat  Transfer 

The  partial  differential  equation  which  describes  the  conduction  of  heat  in  a solid,  when  there  is  no 
radiation  heat  transfer,  does  so  by  giving  the  rate  of  change  of  temperature  at  every  point  within  the  solid. 
It  is  (see,  for  example,  [3]) 
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= divCKgradr)  +5 


(12) 


where  T,  which  is  a function  of  position  and  of  time,  is  the  variable  for  which  we  are  trying  to  find  a 
solution.  The  other  symbols  in  this  equation  are:  p is  the  density  of  the  solid,  c its  specific  heat,  k its 
thermal  conductivity,  and  S is  any  volumetric  heat  source  or  sink,  p,  c,  k,  and  S may  vary  with  position  1 
and  with  T.  Because  these  may  vary  in  Eq(12),  it  is  conceptually  simple  to  consider  a substrate  which  \ 
consists  of  layers  of  material.  Indeed,  it  will  be  straightforward  to  implement  this  in  practice,  as  well. 
Note  that  because  we  only  have  derivatives  of  T,  T need  not  be  the  absolute  temperature:  it  may  be  taken 
to  be  that  relative  to  a convenient  reference  temperature.  In  other  places,  such  as  in  Eq(2),  it  must  be 
taken  to  be  absolute  (i.e.,  in  degrees  Kelvin). 

The  first  term  on  the  right-hand  side  describes  the  diffusion  of  heat  in  the  solid.  If  there  is  any  release 
of  heat  - usually  by  combustion  - at  the  (interior)  point  in  question,  it  is  given  by  S(x,y,z,t).  A general 
expression  for  S is 


(13) 


where  Rj  is  the  reaction  rate  (in  kg/m^-s)  and  is  the  heat  of  combustion  (in  J/kg),  for  the  ith 
reaction.  If  there  is  endothermic  pyrolysis  or  evaporation,  we  have  the  last  term  on  the  right-hand  side 
as  well,  where  is  the  latent  heat  of  evaporation  or  pyrolysis.  An  explicit  expression  for  Rj  will  be 
given  in  Section  2b. 

Since  the  furniture  (apart  from  the  frame)  consists  of  a fabric-covered  pad,  it  is  clear  that  the  program 
must  take  at  least  two  layers  (with  different  properties)  into  account.  Therefore  the  program  was  written 
so  as  to  permit  different  values  for  the  relevant  thermophysical  constants  p,  c,  and  k in  each  layer.  In 
fact,  generally  there  is  not  perfectly  intimate  thermal  contact  between  the  fabric  covering  and  the  padding; 
there  is  a small  but  sometimes  significant  intervening  air  gap.  Normally,  one  would  place  a node  within 
this  gap,  in  order  to  take  a third  layer  into  account;  because  of  the  thinness  of  the  gap,  and  other  technical 
difficulties,  however,  a different  treatment  of  the  effect  of  this  air  gap  has  been  devised:  the  gap  can  be 
represented  in  terms  of  its  "thermal  resistance"  (see  Section  3c.).  This  has  been  programmed  and 
successfully  tested. 

In  writing  Eq.(12),  we  have  made  the  simplifying  assumption  that  the  cushion  is  totally  opaque;  that  is, 
there  is  no  radiative  transfer  of  heat  through  the  cushion.  If  there  were,  the  equation  describing  heat 
transfer  through  the  solid  would  become  still  more  complicated.  However,  the  fabric  and  foam  can  each 
be  thought  of  as  porous,  consisting  of  solid  parts  interspersed  with  void  spaces.  Then  taking  forward 
radiation  transfer  in  those  spaces,  it  is  possible  to  incorporate  a first  approximation  to  (one-d)  radiation 
transfer,  as  shown  by  Kunii  [4] 


K(r)  = + 


(14) 


where  is  the  void  fraction,  and  and  are  the  solid  and  gas-phase  thermal  conductivities, 
respectively;  they  are  each  functions  of  T.  Dp  is  the  mean  pore  diameter,  and 


h,  = 4eaT^ 


(15) 
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Eq.(14)  was  Eq.(5-50)  in  refI2].  Finally,  we  must  note  that  at  a given  temperature,  the  thermal 
conductivity  is  proportional  to  the  density  (also  see  Section  3e): 

K{p)  = (pIp^Kip^)  (16) 


2c.  Time  to  Ignition 


An  experimental  test  for  the  time  to  ignition  is  described  in  Section  4.  We  can  readily  make  an  order-of- 
magnitude  estimate  for  this  time,  based  on  Eq(12):  For  the  very  special  case  that  the  substrate  is  initially 
at  the  uniform  temperature  T^,  that  it  is  homogeneous,  isotropic,  inert,  and  semi-infinite,  that  the 
thermophysical  characteristics  p,  c,  and  k of  the  material  are  independent  of  the  temperature,  and  (finally) 
that  the  problem  is  one-dimensional  (i.e.,  the  incident  flux  is  the  same  everywhere  on  the  plane  z = 0, 
resulting  in  the  slab  being  subjected  to  the  uniform  net  heating  flux  ^netW),  then  Eq(12)  can  be  solved 
analytically  and  explicitly,  and  it  can  be  shown  ([3],  p.76)  that  the  surface  temperature  is  given  by 


T,it)  = + 


^TlKpC 


0 


(17) 


For  the  still  simpler  (unrealistic)  case  = constant,  this  can  be  immediately  integrated,  and  we  obtain 


m = + 


v/UKpC 


(18) 


When  the  left-hand  side  reaches  the  ignition  temperature  Tjg,  the  time  elapsed  must  be  the  time  in  the 
right-hand  side  of  the  equation,  and 


^ _ TlKpC 

A 


/ T -T  \2 
^ o 


4», 


net 


m 


Although  not  a single  one  of  the  simplifying  assumptions  required  to  obtain  Eq(19)  is  valid,  this 
expression  is  nevertheless  very  useful  as  a guide  to  the  form  of  the  dependence  of  t^g.  In  particular,  we 
will  use  Eq(19)  in  order  to  make  sense  out  of  the  experimental  observations,  in  Section  5a. 

(We  can,  if  desired,  get  a somewhat  more  realistic  result  by  relaxing  one  of  the  assumptions  above:  it 
is  possible  to  write,  starting  with  Eq.(17),  an  explicit  expression  for  7^(0  for  the  case  where  the  net  flux 
is  not  constant,  but  results  from  some  impressed  external  flux,  diminished  by  a Newtonian  (i.e.,  linear) 
cooling  loss;  that  is,  where  we  can  write  = h*(Ts  - T^);  see  ref[5]). 


2d.  Pyrolysis  of  the  substrate 

As  described  above,  we  assume  that  the  substrate  consists  of  a thin  fabric  covering  a relatively  thick  foam 
pad,  with  a very  thin  air  gap  between  them.  Each  of  the  two  materials  will  heat  up  and  can  pyrolyze. 
In  this  section  we  describe  how  these  reactions  are  calculated.  Pyrolytic  reactions  can  generally  be 
expressed  in  the  Arrhenius  form 
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kg/m^-s 


(20) 


V^.p/.p.)  =^p7pox^ 


-EJFT 


where  Pf  is  the  density  of  the  fuel  and  p^x  that  of  oxygen;  is  the  activation  energy  for  the  reaction, 
and  R is  the  universal  gas  constant  (not  to  be  confused  wiA  the  reaction  rate  Rp  in  Eq(20)).  It  is 
sometimes  convenient  to  define  an  "activation  temperature"  T^j^: 

= EJR  (21) 


In  the  remainder  of  the  paper,  R is  used  only  to  represent  a reaction  rate.  All  the  quantities  in  Eq.(20) 
except  the  preexponential  factor  A are  at  the  location  (x,y,z)  within  the  porous  solid  and  at  the  time  t; 
this  dependence  is  not  shown  explicitly  in  order  to  minimize  the  complexity  of  the  expression.  Similarly, 
the  subscript  i,  corresponding  to  ith  "reaction"  has  been  suppressed  from  R,  A,  pf,  m,  n,  and  E^.  If  the 
reactions  were  truly  elemental  Arrhenius  reactions,  then  m = n = 1.  However,  Eq.  (20)  is  really  a 
model  equation,  and  therefore  m and  n are  purely  empirical  parameters  (just  as  A and  E^  are),  and  need 
not  be  integers. 

Oxidative  as  well  as  non-oxidative  pyrolysis  (thermal  decomposition)  of  materials  have  been  included. 
The  former  reactions  are  generally  exothermic,  while  the  latter  are  generally  endothermic.  There  may 
be  several  oxidative  reactions  as  well  as  several  non-oxidative  ones;  char  oxidation,  when  it  takes  place, 
is  the  last  step.  In  this  model  a maximum  of  three  reactions  is  permitted  for  each  material.  The  oxidative 
reactions  are  exothermic;  we  assume  that  the  rate  of  each  step  (reaction)  is  adequately  described  by  an 
Arrhenius  equation  of  the  form  (20). 

It  had  been  hoped  that  the  different  pyrolysis  steps  take  place  at  sufficiently  different  temperatures,  that 
the  reaction  sequence  for  the  ignition  of  the  fabric  could  be  taken  to  consist  of  a small  number  of  steps 
which  follow  each  other  sequentially.  That  would  correspond  to  the  equations  taking  the  form 

^ (22) 


where  p^  is  the  density  of  the  fuel  at  the  end  of  the  reaction  step, 
gaseous  species  which  escapes,  the  reaction  rate  is  given  by 

^ dt 

More  generally,  if  species  i is  reacting  and  producing  species  j,  then 


ffi  - ff  - 
dt  ' 3t 


When  the  reaction  produces  a 

(23) 

(24) 


We  will  here  discuss  only  the  pyrolysis  of  the  fabric,  which  is  of  principal  interest:  when  smoldering 
ignition  is  produced,  it  is  almost  always  first  initiated  in  the  fabric.  Pyrolysis  of  the  foam  is  discussed 
briefly  in  Section  4b.  The  fabric  of  particular  interest  is  cotton  duck;  cotton  is  principally  cellulose. 

Analysis  of  the  reaction  kinetics  of  a cellulosic  paper  by  Kashiwagi  and  Nambu  [6]  showed  that  it  could 
be  described  by  three  reactions,  all  of  which  will  proceed  simultaneously  (though  obviously  at  different 
rates).  That  means  that  the  equations  describing  the  creation  and  destruction  of  different  fuel  species  must 


8 


be  included.  Virgin  material  goes  to  char,  via  two  pathways:  degradation  and  oxidative  pyrolysis;  each 
gram  produces  n^  grams  of  char.  The  char  then  goes  to  ash,  via  char  oxidation;  each  gram  of  char 
produces  n^  grams  of  ash.  Thus  one  gram  of  material  produces  n^nj.  grams  of  ash.  The  governing 
equations  are: 

Pv  = (25a) 


and 


The  density  of  the  solids  is 


Pa,  = 

~ o a 


CO 


Hence 


Ps  = Pc*  Pa 


The  subscripts  denote: 
a = ash 
c = char 

CO  = char  oxidation 
d = degradation 


op  = oxidative  pyrolysis 
s = solid 

V = virgin  material 


(25b) 

(25c) 

(25d) 

(25e) 


Besides  the  surface  heat  source  described  in  section  2a,  there  are  also  internal  heat  sources  (and  sinks), 
given  by  Eq(13).  The  heats  of  combustion  must  be  supplied  by  the  user  of  the  model;  those  for  cotton 
are  given  in  sections  4b. 3 and  4b. 4.  The  global  kinetic  constants  given  by  Kashi wagi  and  Nambu  are 
listed  in  Section  4.b.l.  If  the  reaction  rate  of  one  or  more  exothermic  reactions  becomes  high,  for  any 
reason,  the  rate  at  which  energy  is  being  generated  will  exceed  the  rate  at  which  it  is  lost,  and  a "thermal 
runaway"  will  ensue  (very  similar  to  a chain  reaction!).  During  the  runaway,  the  rate  of  pyrolysis  (and 
of  heating)  is  limited  by  the  availability  of  oxygen.  Since  O2  has  0.001  times  the  density  of  the  local 
fuel,  it  would  immediately  be  exhausted,  but  for  the  supply  which  (a)  diffuses  in  from  adjoining  cells, 
and  that  which  (b)  diffuses  in  (or  is  convected  in)  from  the  surface. 


2e.  Gas  Diffusion 

As  oxygen  is  depleted  in  the  regions  where  combustion  (oxidative  pyrolysis,  char  oxidation)  is  taking 
place,  the  concentration  gradient  which  results  will  induce  diffusion  of  oxygen  into  those  regions,  from 
the  surface  as  well  as  from  adjacent,  oxygen-rich  regions,  assuming  the  medium  is  porous  and  permits 
diffusion.  Similarly,  the  gaseous  (and  other)  products  which  are  generated  ~ mainly  CO2  — build  up  in 
local  concentration,  and  diffuse  away. 

The  equations  which  describe  the  rate  of  change  of  species  concentration  are  analogous  to  Eq.  (12)  for 
the  diffusion  of  heat;  assuming  no  convection,  the  ith  species  density  is  governed  by 

= div  (D^  grad  ) + S,  (26) 

01 

where  is  the  diffusion  coefficient  for  species  i in  the  background  "0",  and  Sj  is  the  source/sink  term. 
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If  the  diffusion  coefficient  increases  with  T,  as  it  normally  does,  then  the  CO2  diffusion,  for  example, 
will  take  place  preferentially  towards  the  hot  regions  --  that  is,  towards  the  surface.  The  migration  of 
gases  is  expected  to  have  a minor  effect  on  the  heat  transfer,  and  hence  on  the  temperature  distribution, 
from  this  effect.  However,  the  diffusion  of  oxygen  is  important;  the  rate  at  which  oxygen  enters  the 
reacting  region  sets  the  limit  on  the  reaction  rate. 

A porous  medium  is  one  consisting  of  solid  particles  embedded  in  a gaseous  medium.  Put  in  a different 
way,  a gas  molecule  cannot  penetrate  any  of  the  solid  phase,  but  is  able  to  traverse  the  entire  medium, 
either  because  of  a pressure  difference,  or  merely  from  the  "random  walk"  which  constitutes  diffusion. 
The  diffusivity  of  a gas  species  in  such  a medium  is  really  the  diffusivity  of  those  gas  molecules  through 
a gaseous  medium  consisting  of  the  same  or  some  other  gas,  where  the  solid  volumes  are  excluded.  That 
is  the  effective  diffusivity. 

It  has  been  found  experimentally  (Szekely,  ref[7])  that,  approximately, 

D,„  = Z)„*"  C7) 

where  €>  is  the  fractional  void  space;  this  is  also  referred  to  as  the  "porosity"  in  the  literature.  The 
relationship  is  strongly  dependent  on  the  structure  of  the  (granular)  material;  thus.  Fig.  2.4  from  ref  [7] 
shows  that  for  particles  of  mica, 

O.ffID,  - <27a) 

whereas  for  sand,  a bed  of  glass  spheres,  carborundum  powder,  and  table  salt, 

= 0.677  (27b) 

for  < 0.7.  As  shown  in  Appendix  E,  $ = 0.6  for  the  cotton  fabric  that  will  be  our  main  focus  of 
interest;  hence  Eq.(27b)  is  the  relationship  to  be  used. 

Here  is  the  diffusivity  of  O2  in  N2.  From  Eq.(16.3-1)  of  Bird  et  al  (ref  [8]),  we  find  that  for  O2  in 
N2,  the  temperature  dependence  is 

/ 'p  \1.823 

DTT)  = 0.199  — —\  cmVsec  (28) 

" U93.16j 


2f.  Boundary  Conditions  for  Gases 

There  are  two  cases  of  interest:  first,  when  the  air  above  the  substrate  is  quiescent,  and  second,  when  -- 
as  was  the  case  in  the  experiment  to  be  described  in  Section  4 — there  is  a stream  of  air  impinging  on 
the  substrate.  We  first  obtain  a general  expression; 

When  the  air  above  the  substrate  is  quiescent,  and  there  are  oxidative  reactions  taking  place  in  the 
substrate,  the  concentration  of  oxygen  molecules  in  the  air  decreases  towards  the  surface,  through  a 
boundary  layer.  The  maximum  possible  reaction  rate  is  obtained  by  assuming  y^  = 0,  as  can  be  seen 
from  Eq(29).  On  the  other  hand,  this  is  clearly  not  possible,  since  we  must  have  a finite  concentration 
at  the  surface  in  order  to  get  any  reactions  at  all.  It  is  possible  to  obtain  a value  for  y^  analytically,  if 
one  makes  the  simplifying  assumptions  of  an  isothermal  substrate  and  a first-order  reaction  rate  (ref  [8a]). 
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Mass  transfer  is  a (molecular)  transport  phenomenon,  just  as  is  heat  (momentum)  transfer.  Thus, 
analogous  to  the  first  term  on  the  right-hand  side  of  Eq.(4),  the  mass  transfer  of  oxygen  across  the 
surface  of  the  solid  can  be  written  as 


mo2  ~ 


(29) 


where  is  the  ambient  oxygen  mass  fraction  and  that  at  the  solid  surface;  the  x,y  dependence  of  the 
latter  has  been  put  in  explicitly  in  order  to  emphasize  the  origin  of  the  spatial  dependence  of  m 


The  velocity  at  which  oxygen  enters  the  surface  is 

(see  refl2],  p.  159-160)  where  D is  the  diffusion  coefficient  inside  the  solid.  Thus 


Kas  = P.Y 


In  TMPSUB  (see  [2],  Eq.(5-64),  from  ref  [9])  the  following  equation  was  used: 


Yj,  = 6.38x20 


-3 


7’2.75(j  _ j ^ 123.6) 


RT. 


1/4 


cm/sec 


(30) 


(31) 


(32) 


On  the  other  hand,  in  this  paper  we  obtain  y in  a different  way;  we  will  afterwards  use  Eq(32)  in  order 
to  compare  the  results.  The  way  y is  obtained  is  as  follows:  Because  of  the  similar  origin  of  mass  and 
heat  transfer,  one  can  often  use  the  Reynolds-Colbum  analogy  [10].  It  is  not  difficult  to  show,  from  the 
treatment  in  Chap. 3 of  [10],  that  the  Reynolds-Colburn  analogy  leads  to 


hPr^ 

pc^5c^ 


m/sec 


(33) 


where  h is  the  heat  transfer  coefficient  appropriate  to  the  problem,  Pr  is  the  Prandtl  number.  Sc  is  the 
Schmidt  number  (Sc  = pfD),  and  v is  the  kinematic  viscosity. 


Case  A.  Quiescent  air. 


Here  we  have 


h = kNu/I^  (34) 

where  k is  the  thermal  conductivity,  is  a characteristic  dimension  for  the  problem,  and  Nu  is  the 
Nusselt  number.  We  can  thus  rewrite  Eq.(33)  as 


Y = 


(35) 


In  quiescent  air,  there  is  a constant  movement  of  molecules  in  all  directions;  the  flow  in  any  one  direction 
is  exactly  compensated  by  the  flow  in  the  opposite  direction,  normally.  If,  however,  the  gas  is  near  a 
boundary  which  "absorbs"  some  of  these  molecules  (as  is  the  case  for  oxygen  impinging  on  a reactive 
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substrate)  then  evidently  the  return  flow  is  smaller,  and  hence  the  impinging  flow  is  not  completely 
nullified.  The  result  is  that  there  is  a gradual  decrease  in  the  concentration  of  those  molecules  as  the 
surface  is  approached.  This  region  is  referred  to  as  the  "boundary  layer." 

The  quantities  in  Eq.(35)  are  to  be  evaluated  at  some  characteristic  (mean)  temperature  of  the  boundary 
layer.  That  temperature  can  be  taken  to  be  the  mean  between  T^  and  T^.  For  quiescent  air,  with  the 
surface  horizontal  and  facing  upward,  /<,  is  the  boundary  layer  thickness.  The  expression  for  the  Nusselt 
number  is  given  in  section  2a. 

Case  B.  Impinging  air 

The  Nusselt  number  for  this  case  is  given  by  Eq.(D10),  Appendix  D.  Inserting  that  into  Eq.(35),  we 
obtain 

Y - 0.767  Sc 

where  a is  the  thermal  diffusivity  of  air.  Re  = vlJJv  is  the  Reynolds  number  and  u^,  is  the  incoming 
velocity  of  the  impinging  flow. 

We  now  show  that  when  oxygen  consumption  is  significant,  it  takes  place  in  a very  narrow  layer;  this 
is  analogous  to  the  flame  sheet  approximation  for  combustion  in  air.  As  a result,  we  may  take  the 
concentration  of  oxygen  in  that  layer  to  be  some  appropriate  average  value,  rather  than  having  to  actually 
program  in  the  species  diffusion  equations,  via  Eq(26). 


2g.  An  Approximation  to  Oxygen  Diffusion 

In  this  section,  it  is  shown  that  a reasonable  approximation  to  the  overall  reaction  rate  can  be  obtained 
without  explicitly  solving  the  gas  diffusion  equation,  Eq.(26),  for  the  diffusion  of  oxygen  in  the  substrate. 
This  is  the  approximation  made  in  the  model. 

In  the  experiment  to  be  discussed  in  Section  4,  a jet  of  air  is  directed  downward  at  the  substrate.  Hence 
there  is  no  boundary  layer  in  this  case,  and  the  oxygen  concentration  just  outside  the  surface  is  21%. 
When  the  reaction  rate  is  low,  the  oxygen  concentration  in  the  top  of  the  substrate  will  not  be  much 
affected  by  the  slow  reactions.  When  the  reaction  rate  becomes  high,  on  the  other  hand,  the  mean 
oxygen  concentration  in  that  thin  layer  becomes  low;  it  cannot  get  too  low,  however,  since  (as  discussed 
in  Section  2f),  the  reaction  rate  would  then  fall.  Hence  a quasi-steady  state  is  established.  Thus,  in  the 
"runaway"  phase,  the  reaction  rate  must  be  given  by  the  boundary  condition  Eq(29). 

We  will  now  argue  that  when  significant  oxidative  pyrolysis  is  taking  place, 

(a)  taking  a constant  value  for  the  oxygen  concentration  ([O2])  in  the  substrate  is  a valid  approximation, 

(b)  most  of  the  pyrolysis  occurs  in  the  surface  layer,  and  that  therefore 

(c)  we  obtain  an  adequate  approximation  to  the  pyrolysis  rate,  without  explicitly  including  the  oxygen 
diffusion  equations: 

When  the  surface  temperature  is  relatively  low,  there  is  negligible  pyrolysis,  and  the  oxygen  distribution 
[02](x)  in  the  substrate  is  approximately  uniform.  When  Tg  has  risen  to  the  point  that  perceptible 
pyrolysis  is  taking  place  at  and  near  the  surface,  the  [O2]  profile  in  the  substrate  dips  as  the  surface  is 
approached  from  below  (i.e.,  from  within);it  must  reach  a minimum  and  then  rise  again  at  the  surface. 
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because  of  the  diffusion  of  oxygen  from  outside.  When  the  temperature  reaches  values  so  high  that  [OJ 
is  pulled  down  to  negligible  values  near  (but  not  the  surface,  the  principal  source  of  oxygen  in  the 
reaction  zone  is  from  the  air  diffusing  in  from  the  surface;  the  region  where  there  is  significant  oxidative 
reaction  is  then  highly  localized  near  the  surface,  with  [O2]  being  highest  at  the  surface,  and  falling 
rapidly  until  it  reaches  negligible  values.  The  characteristic  oxygen  penetration  distance  is  6,  and  it  is 
clear  that  6 is  a steeply  falling  function  of  T.  (At  some  further  distance  in,  the  oxygen  must  diffuse 
towards  the  surface  from  other  (deeper,  or  peripheral)  parts  of  the  porous  substrate,  because  of  the 
concentration  gradient,  so  that  the  profile  must  rise  again  as  we  go  deeper  --  it  can  do  so  because  at  that 
depth  the  temperature  has  fallen  sufficiently  to  "freeze  out"  the  (oxidative)  reaction  rate  --  i.e.,  it  is 
negligibly  small). 

Suppose  that  is  high  enough  to  drive  [O2]  to  near-zero  at  some  depth.  Then  we  may  write,  crudely, 

Yix)  - (37) 

with  Yg  = surface  value  of  [O2],  and  6 = characteristic  penetration  depth.  If  we  write  the  reaction  rate 
in  the  form 


R[x,nx)]  = y(jc)"F[jc,7’(x)]  = y"G(x)  (38) 

then  the  total  reaction  rate  per  unit  area  is 

w"  = fy(xyGix)dx  (39) 

J o 

where  G(x)  is  a fairly  steep  (decreasing)  function  of  x.  We  may  further  write  (still  very  crudely) 

G{x)  - (40) 

Then 


We  can  also  use  the  mean-value  theorem,  and  write  the  integral  in  Eq(41)  as 


Then 


rY"G^e-^'f^dx  = Y”G^Qs/^l2 

J o 


n0 


,2\ 


-1/2 


(41) 


(42) 

(43) 


We  now  need  an  estimate  for  6/6.  Clearly,  the  reaction  rate  must  be  significant  for  a distance  comparable 
to  the  penetration  depth  6:  else,  if  the  reactions  "froze"  before  [O2]  becomes  negligible,  [O2]  would  not 
fall  below  some  appreciable  value.  Moreover,  the  reaction  rate  cannot  be  significant  below  that  depth, 
precisely  because  there  is  so  little  oxygen  below.  Thus  6 = 8.  For  n,  we  take  n,,  = 0.5,  as  given  in 
[6].  Then  Eq.(43)  yields 

Y = 2YJ3  (44) 

Next,  consider  the  effect  of  discretizing  the  equations  (see  section  3).  If  6 > Ax,  then  assuming  that  the 
reaction  takes  place  in  the  top  cells  only,  clearly  underestimates  the  reaction  rate.  However,  when  8 is 
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large,  the  total  reaction  rate  is  very  low.  That  is,  while  the  fractional  error  is  large,  the  absolute  error 
is  small,  and  (we  shall  show)  typically  negligible. 

Assume,  for  the  sake  of  simplicity,  only  one  oxidation  reaction.  Assume,  also,  that  the  temperature  and 
reaction  rate  are  uniform  in  a thin  sheet  of  depth  6 at  the  surface,  where  T is  maximum  (this  is  essentially 
what  is  implicitly  assumed  in  the  numerical  calculation).  Then  if  the  stoichiometric  fuel/oxygen  ratio  is 
r,  the  fuel  bumup  rate  is 

R6  = m'l  = (45) 


where  R is  the  reaction  rate  (g/cm^sec  or  kg/m^s).  If  6 is  greater  than  the  thickness  of  the  surface  cell, 
then  evidently  the  reaction  rate  in  the  top  layer  is  limited  to  RAz;  whereas  if  6 < Az,  Eq.(45)  gives  the 
limit.  Thus  the  reaction  rate  in  the  top  cell  is 


R^opAz  = min[RAz,  rm^"]  g/cm^-sec 


(46) 


Note  that  when  the  reaction  rate  is  so  high  that  6 < Az,  the  concentration  at  the  bottom  face  of  the  top 
cell(s)  is  close  to  zero;  therefore  a reasonable  value  to  take  for  [O2]  in  the  cell  is  on  the  order  of  the  mean 
value,  11%,  (or  0.12  for  the  mass  fraction);  or,  according  to  Eq.(44),  < Y>  = %(0.2318)  = 0.15 

We  now  make  a calculation  to  find  a typical  value  for  6.  As  is  clear  from  Eq.(45),  we  must  begin  by 
calculating  a reaction  rate  R.  In  order  to  find  R,  we  must  anticipate  some  results  from  section  4b. 3: 
using  Eqs(79)  and  (80),  we  find  that 

Rop  = Pvkop  = Pv  1.5x10'“  Y« " (WjAV„)'-3  exp(-160/RT) 
and 

R=o  = pA.  = Pc  3.4x10"  (W.AVJ  exp{-160/RT) 

The  meaning  of  the  subscripts  is  made  clear  in  the  list  below  Eq.(25e).  These  rates  are  in  kg/m^-min. 
We  must  choose  characteristic  values  for  Py,  Y,  W^j/W^,  p^.,  W^.AVq,  and  T,  in  order  to  get  estimates. 
Reasonable  values  are 

Wj  = Wo/2,  Wo  = Wo/4,  and  Y = 0.11 

For  Py  and  p^,  we  assume  that  the  volumes  are  unchanged,  so  that 

Py/Po  = WdAVo  and  = 

Finally,  we  take  p^  = 1560  kg/m^.  We  then  find  that 

kop  = 2.02x10'^  exp(-160/RT)  min'* 
and 

koo  = 1.52x10*°  exp(-160/RT)  min'* 

The  universal  gas  constant  R = 8.31447  J/mol-K,  so  that  the  activation  temperature  for  both  reactions 
is  = 160,000/R  = 19244  K.  For  a temperature  in  the  vicinity  of  the  ignition  temperature,  say, 

400°C  = 673.16  K,  we  then  have 
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yielding 

and 


op 


7.77  min* 


and  kgo  = 0.00584  min’*, 


I 

! 


R 

R 


op 

CO 


/^V^op 

“ Pc^co 


101  kg/m^-s 
0.038  kg/m^-s. 


Thus  at  400 °C,  with  these  assumptions  for  and  p^,  oxidative  pyrolysis  takes  place  more  than  2600 
times  faster  than  does  char  oxidation,  and  we  may  therefore  take 

R = Rop  = 101  kg/m^-s 


Next,  we  need  to  have  the  rate  at  which  oxygen  enters  the  medium,  . Eqs(29)  and  (31)  give  the 

required  value.  For  the  ignition  experiment,  the  appropriate  expression  to  use  to  find  y is  Eq.(36).  The 
temperature  at  which  the  various  quantities  must  be  evaluated  is  that  of  the  purging  gas.  Examination 
of  Table  D1  in  Appendix  D shows  that  Tg  is  quite  high.  Assume  Tg  = 800  K.  Eq(28)  yields  0^(800) 
= 1.24  cm^/sec.  Tlien  the  Schmidt  number  is  Sc  = 0.663,  and 

7 * 8.10x10’^  vfe  m/s 

(with  fj.  in  meters).  The  characteristic  air  velocity  for  calculating  Re  is  given  by  Eq(D15)  in  Appendix 
D;  finally,  the  characteristic  length  is  the  standoff  distance,  6.  With  6 = 5.4  mm,  we  find  Re  = 1 1 .06 
and  7 = 4.99  cm/sec. 

(Incidentally,  if  we  had  quiescent  air  at  800  K,  Muaramatsu’s  expression,  Eq(32),  yields  7jj  = 2.35 
cm/sec). 

Finally,  taking  = 0.232  and  = 0,  Eqs(29)  and  (31)  then  yield 

« 13.6g/m2-s 

We  must  also  have  r (the  stoichiometric  fiiel/oxygen  ratio),  in  order  to  use  Eq(45)  to  get  the  penetration 
depth  6.  The  stoichiometric  air/fuel  mass  ratio  for  wood  is  S = 5.78,  close  enough  to  that  for  cellulose. 
Hence  r = 0.746,  and  we  find  6 = 0.10  mm.  Thus  the  penetration  depth  b is  indeed  much  smaller  than 
the  layer  thickness  Ax  (which  is  of  the  order  0.5  mm). 

For  Tg  = 300°C,  on  the  other  hand,  we  find  k^p  = 0.053  min’*,  so  that  R(,p  = 1.38  kg/m^-s.  Then  6 
« 7.3  mm,  several  times  Ax,  but  the  resulting  mass-loss  rate  is  only  0.043  g/m^-s,  which  is  indeed 
negligible. 

For  the  case  of  quiescent  air  above  the  substrate,  the  mean  value  for  the  oxygen  concentration  to  be  used 
in  the  top  cell  during  the  runaway,  is  midway  between  the  value  at  the  top  surface  (i.e.,  at  the  bottom 
of  the  boundary  layer),  and  at  the  bottom  of  the  top  cell  (presumably,  close  to  zero). 

The  only  time  that  we  cannot  justifiably  make  the  simplification  that  the  oxygen  concentration  in  the  top 
layer  is  either  that  at  the  top  surface  or  half  (or  2/3ds)  that  value  (during  runaway),  is  when  the  reaction 
rate  is  intermediate  between  the  very  low  values  and  the  runaway  value.  This  period  should  be  relatively 
short,  and  the  error  introduced  by  these  simplifications  should  not  be  large. 
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3.  NUMERICS 


In  general,  it  is  not  possible  to  solve  Eq.(12)  analytically,  in  spite  of  the  simplifying  assumption  of  no 
radiation  heat  transfer,  and  even  for  the  case  S=0.  That  is,  to  write  down  an  explicit  expression  which 
gives  T in  terms  of  the  inputs.  This  is  so,  because  of  the  nonlinear  and  nonuniform  boundary  conditions. 
It  is  therefore  necessary  to  resort  to  a numerical  procedure,  which  is  that  incorporated  in  TMPSUB2. 

3a.  Introduction 

The  development  of  TMPSUB2  centers  around 
the  capability  to  simulate  transient  heat  transfer. 

A one-dimensional  heat  conduction  problem 
provides  the  simplest  example  to  illustrate 
transient  simulation  methods.  Figure  1 shows  a 
portion  of  a one-dimensional  conduction  problem 
in  which  the  material  has  been  divided  into  thin 
layers.  This  example  will  be  described  by 
physical  instead  of  mathematical  arguments 
following  the  description  given  by  Clausing 
(refill],  pp.  157-213). 

The  figure  focuses  on  a representative  layer  of 
material  of  thickness  Ax;  centered  in  a node  at 
coordinate  Xj.  This  material  layer  is  represented 
by  a single  temperature  T;  and  a corresponding 
thermal  conductivity  Xj,  density  pj,  and  specific 
heat  Cj.  Assume  this  layer  has  a surface  area  of  magnitude  A in  the  Y-Z  plane.  The  distance  between 
nodes  i and  i-1  is  given  by  (=  Xj  - Xj.j).  Subscripts  i and  n refer  to  positions  in  space  and  time, 
respectively.  Subscript  n is  temporarily  suppressed  until  it  becomes  necessary  to  consider  time,  in  the 
following  equations. 

The  instantaneous  internal  energy  of  layer  i is  given  by: 

U.  = pc  A FT.  = pc  Ax,.  AT.  = C.J,.  (47) 

where  T is  the  absolute  temperature.  This  equation  also  serves  to  define  the  heat  capacity,  C,  assigned 
to  node  i.  Heat  is  transferred  to  and  from  layer  i by  three  methods: 

(1)  conduction  from  layer  i-1  at  the  rate 

9,.  = ( T;.,  - r,)  k,.  X / 6a:,  = K,_  ( T,.,  - T,)  (48) 

(2)  conduction  from  layer  i 4- 1 at  the  rate 

= ( Ti.i  - 7’i)  K,.  A I - T.)  (49) 

and  (3)  internal  heat  generation  or  radiation  absorbed  within  the  layer,  Sj. 

Eqs(48)  and  (49)  define  the  thermal  conductance,  K,  in  each  direction.  The  possibility  of  a thermal 
conductivity  which  varies  with  position  is  explicitly  taken  into  account  here  using  labels  -I-  and  — : 
and  xj+.  A good  first  order  approximation  for  K■^_  is  the  harmonic  average  of  x;  and  Xj_j: 


layer  I 

AX|  — 

^ ^ Y 

1-^—6: 

^i-1 

ax  i I 

2 

<i — 

• A 

Figure  1 . One  Dimensional  Conduction 
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(50) 


2K.K..1 


Ki+K..i 


with  a similar  expression  for 


The  change  in  internal  energy  of  this  layer  between  time  t„  and  time  t^^.!  is  given  to  first  order  in  At  by 

^^,.-..1  " = Arc, (51) 


where  the  time  dependence  has  now  been  put  in  explicitly.  There  are  several  common  solutions  to  Eq. 
(51)  depending  on  when  the  heat  gains  are  evaluated.  One  solution  involves  evaluating  at  time  n: 


C T 


(52) 


which  is  the  standard  Euler  explicit  time  integration  formula.  "Explicit"  means  that  T|  ^,4.)  can  be 
directly  computed  from  values  known  at  time  n.  On  the  other  hand,  evaluating  at  time  n + 1 gives; 


C T 


= ^ AT,., 


k,AT, 


»+'•  i+l,«+l 


-T 


i.n+l 


(53) 


which  is  Euler’s  standard  implicit  time  integration  formula.  "Implicit"  means  that  Tj^+j  is  computed 
from  other  values  also  evaluated  at  time  n+ 1.  These  values  depend  implicitly  on  each  other  and  must 
be  computed  by  a solution  of  simultaneous  equations. 


Clausing  ([11],  p.  190)  also  gives  a discussion  of  stability  in  terms  of  thermodynamic  laws.  Rearranging 
Eq.(52)  to  solve  for  Tj  gives 


K.  r. , +5.  )Ar/c.  ,+r.  [C.  -(K. 

I+,B  l + l,B  (.B'  ' «,B+1  I.B*^  l,B  l-.r 


(54) 


where  the  time  subscript,  n,  has  been  added  to  the  K terms  to  indicate  exactly  when  these  values  are 
evaluated.  For  the  sake  of  argument,  assume  = Cin(=Cj)for  simplicity.  There  is  no  solution 
if  Cj  = 0.  If  C,  is  sufficiently  small  or  At  sufficiently  large,  then  ( Ki_  „ + Kj^.  „ ) At  / Cj  > 1,  and 
as  Tj  n increases  T;  must  decrease,  and  vice  versa.  This  is  thermodynamically  impossible.  It  shows 
up  in  a numerical  solution  as  oscillations,  i.e.  "instability",  in  the  node  temperatures  at  each  time  step. 
These  oscillations  tend  to  quickly  increase  to  totally  meaningless  values.  In  general,  the  smaller  the 
thermal  mass  of  the  element,  the  smaller  the  time  step  needed  for  a stable  explicit  solution.  This  suggests 
a simple  technique  to  determine  the  minimum  stable  time  step  for  any  element  in  the  system. 


Thus,  rearranging  the  implicit  Eq(53)  to  solve  for  Tj  gives 


i,B+l 


^i,nTi,n  ^ ^ (^<-.B4l  ^-l.B*!  '*~^i*l,B4l  '^'^t.B+l) 


(55) 


This  equation  shows  none  of  the  computational  or  thermodynamic  problems  of  Eq.(54)  indicating  that 
the  standard  implicit  method  is  stable  for  all  time  steps. 

The  spatial  discretization  error  (for  a uniform  grid)  for  the  standard  explicit  and  standard  implicit  methods 
is  proportional  to  (Ax)^  (for  a variable  grid,  the  accuracy  is  reduced  somewhat;  see  Section  3d).  The 
time  discretization  error  is  proportional  to  At. 
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The  standard  explicit  and  standard  implicit  methods  err  in  opposite  directions.  Therefore,  a more 
accurate  solution  can  be  obtained  by  combining  the  two  methods.  Expressing  this  combination  generally 
in  terms  of  a parameter  /3  gives: 


+5.  ) + 1 1+5  ,)] 

^i,n'  + l Ti+.n  + l + 


(56) 


where  0 < /3  < 1. 

/3  = 0 corresponds  to  the  standard  explicit  method, 

/S  = 1/2  corresponds  to  the  Crank-Nicholson  method, 
/S  = 2/3  corresponds  to  the  Galerkin  method,  and 
/S  = 1 corresponds  to  the  standard  implicit  method. 


For  /S  ^ 1/2  this  method  is  unconditionally  stable,  although  the  solution  may  be  oscillatory.  For  j(3  > 
3/4  (approximately)  the  solution  is  stable  and  non-oscillatory.  For  jS  = 1/2,  the  time  discretization  error 
is  proportional  to  (At)^. 


The  methods  presented  above  extend  directly  into  three  dimensions.  For  a Cartesian  coordinate  system 
and  a cell  of  dimensions  Ax;  by  Ayj  by  Az^,  Eq.(51)  can  be  rewritten  to  account  for  conduction  from 
the  six  adjacent  cells  in  the  3-D  system: 


(57) 


where 


WvAy,  Az^/6x. 

qj.  = -r.^.t)V  Ax.  Az^/dy. 

% = (^ij*U-^/j.i)VAx,.Az*/6y^,j 

<lk-  = (’^ij.k-x  - T'ij.k')  ^k-  A^i  Ay.  / hz, 
^k*  = V Ax.  Ay.  / 6zt,i 


and  Sj  j ,5^  represents  other  heat  added  directly  to  the  cell. 


3b.  Boundary  Conditions 

Special  treatment  is  required  for  cells  on  the  boundaries  of  the  region  being  modeled.  In  particular,  the 
nodes  which  represent  the  cells  are  placed  on  the  boundary,  rather  than  at  the  center.  Referring  to  the 
one-dimensional  example  in  Figure  1,  the  surface  layer  is  only  half  as  thick  as  the  others. 

An  adiabatic  boundary  condition  (b.c.)  is  handled  by  setting  the  appropriate  heat  flux  terms  in  Eq(45)  to 
zero.  A constant  temperature  (isothermal)  b.c.  is  handled  by  leaving  the  temperature  unchanged. 

The  surface  of  the  substrate  (z=0  plane;  see  Fig.  3)  transfers  heat  to  the  environment  by  convection  and 
by  radiation.  The  convective  heat  gain  for  cell  i,j,l  is  given  by 
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(58) 


where  is  the  temperature  of  the  surrounding  (ambient)  air,  and  h is  the  heat  transfer  coefficient. 

A positive  value  of  S^.  represents  a heat  gain  by  the  cell.  The  radiative  heat  gain  is  given  by 

where 

is  the  temperature  of  the  surrounding  surfaces, 
a is  the  Stefan-Boltzmann  constant,  and 

e is  the  emissivity  of  the  fabric  (assuming  that  a,  the  absorptivity,  = e). 

The  temperatures  in  Eq(59)  must  be  absolute  (Kelvin)  temperatures.  In  TMPSUB2  the  temperature  of 
the  surrounding  surfaces  is  assumed  equal  to  the  air  temperature.  (Note  that  Eq(59)  can  be  rewritten  in 
an  apparently  linear  form  similar  to  Eq(58); 


(59a) 


This  form  will  be  useful  further  on. 


The  heat  flux  from  a smoldering  cigarette  to  the  fabric  is  represented  by  the  following  pair  of  equations: 

.2 


X*  / 


ix^x„+vt) 


(60a) 


and 


x-x  -vfV 


X- 


y J 


(x<x  +Vf) 


(60b) 


This  heat  flux  is  converted  to  a heat  gain,  S^,  by  multiplying  the  flux  at  the  position  of  the  node  by  the 
cell  surface  area.  These  values  are  adjusted  slightly  so  that  their  (discrete)  sum  is  equal  to  the  total  heat 
flux  represented  by  Eq.(40a,b): 


I'jy.dxdy 


+ o. 


(60c) 


3c.  Air  Gap 

Since  the  furniture  (apart  from  the  frame)  consists  of  fabric-covered  padding,  it  is  clear  that  the  program 
must  take  at  least  two  layers  (with  different  properties)  into  account.  Therefore  the  program  was  written 
so  as  to  permit  different  values  for  the  relevant  thermophysical  constants  p,  c,  and  k at  each  node.  In 
fact,  generally  there  is  not  perfectly  intimate  thermal  contact  between  the  fabric  covering  and  the  padding: 
there  is  a small  but  sometimes  significant  intervening  air  gap.  Normally,  one  would  place  a node  within 
this  gap,  in  order  to  take  a third  layer  into  account;  because  of  the  thinness  of  the  gap,  and  other  technical 
difficulties,  however,  a different  treatment  of  the  effect  of  this  air  gap  has  been  devised:  the  gap  can  be 
represented  in  terms  of  its  "thermal  resistance". 
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The  optional  air  gap  between  the  fabric  and  the 
padding  is  modeled  by  assuming  one-dimensional 
heat  transfer.  Figure  2 shows  the  basic 
configuration  and  nomenclature  for  an  air  gap  of 
thickness  s between  cells  ij,k  and  ij,k+l.  s is 
small  relative  to  the  separation  between  the  cells, 
6^+1.  Although  pyrolysis  of  the  fabric  and  the 
possible  melting  and/or  pyrolysis  of  the  padding 
may  well  change  the  dimensions  of  the  air  gap, 
the  simplifying  assumption  is  nevertheless  made 
here  that  the  air  gap  has  fixed  and  uniform 
dimensions.  The  heat  transfer  between  these  cells 
is  given  by 


Figure  2.  Air  Gap  Heat  Transfer  Model 

(61) 


where  K is  an  implicit  function  of  qij,k+’ 

The  overall  heat  transfer  coefficient  is  given  by 


K = 


1 1 

— + — + ■ 


where  the  fabric  conductance  is 

^/ = 2Ax.Ay.K.^.^/6;t*i  > 

the  padding  conductance  is 

Kp  = 2Ax,A)^k,.^  ^^j/ , 

the  fabric  (bottom)  surface  temperature  is 

~ ^ij.k  ~ ^ij.k*  1 ^u.k  ’ 

the  padding  surface  temperature  is 

the  radiant  conductance  is 

_ q Ayj(T^+Tp)iTf^f^) 


and  the  convective  conductance  is  = Ax.Ay^.  k^/5  . 


The  conductance  of  the  air,  is  evaluated  at  the  average  of  Tf  and  Tp.  These  equations  are  solved  to 
give  K and  qij,k+  during  the  overall  process  used  to  compute  cell  temperatures. 

3d.  Variable  Grid 

The  heat  from  the  cigarette  spreads  into  the  substrate  by  conduction.  In  order  to  correctly  estimate  the 
temperatures  near  the  peak,  it  is  important  that  conduction  to  the  outer  boundaries  (at  x=0,  x=Xn,^, 
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y=y^ax’  negligible,  since  it  cannot  be  known  a priori.  This  requires  a relatively  large 

region  relative  to  the  size  of  the  cigarette  heat  flux  pattern.  However,  setting  Ax,  Ay,  and  Az  sufficiently 
small  to  achieve  the  desired  accuracy  in  the  conduction  calculation  can  result  in  an  extremely  large 
number  of  cells  ( N = (x^j^x/^)  ‘ (ymax^^y)  ‘ ) with  correspondingly  large  memory 

requirements  and  computation  time. 

The  heating  flux  from  a smoldering  cigarette  rises  from  negligible  values  to  a high  peak,  on  the  order 
of  60  kW/m^  over  a region  only  a few  millimeters  in  extent.  In  order  to  follow  this  faithfully,  the  region 
must  be  covered  by  a mesh  which  is  fine  enough  so  that  there  are  no  changes  from  one  mesh  point  to 
another  large  enough  to  produce  numerical  inac- 
curacies or  instabilities.  Thus,  the  required  size 
of  the  grid  is  inversely  proportional  to  the 
temperature  gradient,  where  steep  temperature 
gradients  occur  only  near  the  point  of  peak  heat 
flux.  Therefore,  a variable  grid  is  used.  This 
grid  consists  of  a few  constant-  width  cells  near 
the  peak  followed  by  cells  of  regularly  increasing 
size  to  the  outer  boundaries.  The  increase  in  cell 
size  is  based  on  geometric  progression,  and  can 
be  different  for  each  axis.  Thus,  for  example, 

Ayj+i  = Ry  Ayj  (Ry  ^ 1).  The  general  effect  of 
this  variable  grid  is  shown  in  figure  3.  This  is 
easy  to  implement  in  Eq.(45)  which  explicitly 
incorporates  the  grid  sizes.  The  variable  grid 
gives  results  which  are  less  accurate  than  a 
constant  grid.  A benchmark  test  indicates  that  the 
results  for  R = 1.23  differ  by  < 0.11%  at  t=  100  sec.  from  those  obtained  for  the  constant  grid  case, 
with  still  smaller  errors  for  R closer  to  unity  (see  Table  Al,  Appendix  A). 

Note:  In  choosing  a grid,  the  user  determines  four  items,  in  each  coordinate  direction:  the  width  of  the 
constant-width  cells  in  the  fine-grid  region,  the  number  of  such  cells,  the  total  number  of  nodes  along  that 
axis,  and  the  total  width  in  that  direction.  The  program  then  does  the  arithmetic,  and  finds  the 
corresponding  value  of  Rj.  If  that  is  < 1,  an  error  message  will  be  returned.  Likewise,  if  too  many  grid 
points  result  from  the  attempted  selection,  an  error  message  will  result. 

The  point  of  peak  heat  flux  moves  as  the  cigarette  smolders.  TMPSUB2  adjusts  the  x-coordinates  to  keep 
the  fine  grid  region  centered  on  that  peak.  This  adjustment  is  made  by: 

(1)  computing  the  new  x-coordinates  of  the  shifted  grid, 

(2)  computing  a cubic  spline  curve  for  the  temperature  in  each  row  of  cells  in  the  old  grid,  and 

(3)  using  the  curve  to  compute  the  cell  temperatures  in  each  row  at  the  new  grid  positions. 

The  standard  explicit  algorithm  for  solving  Eq.(12),  without  pyrolysis  included,  was  checked  in  several 
ways,  principally  by  comparing  its  predictions  against  known  analytical  solutions  (see  Appendix  A). 
These  checks  showed  that  the  numerical  procedure,  and  the  computer  program  for  implementing  it,  are 
correct,  effective,  and  accurate. 
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3e.  Variable  Thermal  Properties 


The  conductivity  and  specific  heat  of  the  substrate  materials  are  known  to  vary  with  temperature. 
TMPSUB2  allows  the  user  to  describe  this  variation.  The  values  entered  for  ac  and  c are  converted  to 
cubic  spline  curves  giving  k(J)  and  c(T).  The  conductivity  is  further  adjusted  within  the  program  to  take 
into  account  the  fact  that  the  thermal  conductivity  is  proportional  to  the  density: 


(62) 


K = p k(J)  / Po, 


where  p^  is  the  original  density  (recall  that  the  assumption  of  constant  k,  p,  and  c,  in  section  2c,  was  a 
special  case  used  only  in  order  to  obtain  Eqs(17)  to  (19)). 

3f.  Pyrolysis 

Consider  pyrolysis  from  the  finite  difference  viewpoint.  We  begin  with  a cell  having  constant  volume, 
AV  (=  Ax  Ay  Az)  and  temporarily  ignore  the  subscripts  i,J,k.  This  cell  contains  a mass  Py  AV  of  virgin 
(unpyrolyzed)  material,  a mass  p^  AV  of  char,  and  a mass  p^  AV  of  ash.*  The  numerical  method  must 
keep  track  of  the  density  of  each  material  as  a function  of  time;  in  this  program,  this  is  done  via  the 
difference  equation 


(63) 


Px.n+1  ~ Pjt,B  ^ 9jt,n+P 


where  the  subscript  x may  be  v,  c,  or  a,  and  jS  indicates  the  type  of  time  integration,  as  discussed  above, 
is  modeled  using  Eqs.(25a,b,c).  Note  that  p^  = p^  = 0 and  p^  = p^  = p^  at  t = 0. 

Some  material  is  converted  to  gases  during  pyrolysis;  the  rate  as  which  gases  are  created  is  given  by 
Eq.(25e).  The  gases  are  lost  from  the  cell,  and  they  carry  away  all  the  energy  they  contain  --  i.e.,  there 
is  an  enthalpy  loss  which  does  not  affect  the  temperature  of  the  remaining  mass.  Moreover,  since  the 
gas  diffusion  equations  have  not  been  explicitly  included,  any  possible  loss  of  heat  from  the  escaping  hot 
gases  to  the  cooler  solid  in  other  regions  of  the  substrate  is  ignored.  Thus,  the  rate  of  heat  gain  in  the 
cell  due  to  pyrolysis  is  (almost  exactly)  given  by 


(64) 


The  Rj  are  given  in  Section  2c,  and  dp^/dt  is  given  by  Eq(25e). 

3g.  Time  Integration 

A choice  must  be  made  for  the  time  integration  method.  Following  the  discussion  of  Belytschko  (ref[12], 
pp.55,  419,  445),  the  advantages  of  explicit  time  integration  are; 

(1)  Fewer  calculations  per  time  step. 

(2)  Algorithm  logic  and  structure  are  simple;  this  implies  that  it  is  good  for  testing  new  ideas. 

(3)  Complex  nonlinearities  are  easily  handled. 

(4)  It  requires  little  core  storage  compared  to  implicit  methods  using  direct  elimination  procedures. 

* That  is,  the  densities  are  here  defined  on  a bulk  basis. 
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(5)  It  is  very  reliable  in  terms  of  accuracy  and  completing  the  computation. 

The  only  notable  disadvantage  is  that  explicit  time  integration  is  only  conditionally  stable  so  that  a very 
large  number  of  time  steps  may  be  required. 

With  regard  to  accuracy,  since  implicit  methods  are  unconditionally  stable,  they  can  easily  be  used  with 
too  large  a time  step,  leading  to  significant  time  integration  errors.  The  stability  requirements  for  explicit 
time  integration  force  the  time  step  to  be  so  small  that  the  time  integration  error  is  almost  always  smaller 
than  the  spatial  discretization  error.  Of  course,  it  is  also  possible  to  use  a spatial  discretization  that  is 
much  too  large. 

The  addition  of  pyrolysis  to  the  model  required  significantly  smaller  cells  in  the  region  of  interest,  and 
therefore  required  significantly  longer  execution  time.  Hence  a better  method  than  explicit  time 
integration  was  required.  The  following  method  attributed  to  Saul’yev  as  described  by  Larkin  [13]  and 
Clausing  [11]  was  adopted: 


Again  consider  the  one  dimensional  presentation  of  Fig  1.  Assume  that  the  calculation  of  cell 
temperatures  is  proceeding  in  the  positive  x direction.  Then  at  cell  i,  Tj.j  is  a known  quantity  and 
can  be  used  in  computing  T;  „+j.  Eq.(51)  becomes 


Af/. 


i,n-«+l 


(65a) 


During  the  next  time  step,  calculate  cell  temperatures  in  the  negative  x direction.  In  that  case  Eq.(51) 
becomes 


AU.  , = C.  ,r.  ,-C.  T = AHq. 

i.n+l  i,n*l  i,n  i.n  '’i-,/ 


?<+,«+!  ^i,n+p  ^ 


(65b) 


Both  Eqs(65a)  and  (65b)  are  unconditionally  stable  because  of  the  inclusion  of  the  implicit  terms; 
operating  together,  the  truncation  errors  during  successive  time  steps  tend  to  cancel  leading  to  an  O(At^) 
algorithm.  This  method  is  directly  expanded  into  2 or  3 dimensions  by  adding  the  j and  k position 
indices  and  the  qj.,  qj+,  q^.,  and  heat  gains.  The  key  factor  is  that  it  is  only  necessary  to  solve 
implicitly  for  one  cell  temperature  at  a time.  There  are  no  time  consuming  simultaneous  equations  to  be 
solved.  Tests  indicate  that  this  method  is  very  accurate  except  for  the  possibility  of  some  small 
oscillations  as  with  the  Crank-Nicholson  method. 


A question  remains  on  when  to  evaluate  Sj,  Cj,  and  Xj  (which  is  implicit  in  the  q’s)  in  Eqs(65a). 
Numerical  errors  are  minimized  by  evaluating  S;  at  time  step  n+*/^.  These  terms  are  all  functions  of 
temperature  in  the  TMPSUB2  program.  Furthermore,  they  are  such  complicated  functions  that  the 
equations  cannot  be  solved  directly.  The  following  equation  must  be  solved  implicitly  for  Tj  j „+i. 


^ij,k+,n  ^ 


where 


(66) 


^i.y, t,  m ~ P ij,  k,m  ^ ^i,J,  k k,  m k, m 

with  m = n and  n+ 1.  The  total  solid  density  in  the  cell  is  the  sum  of  the  virgin,  char,  and  ash  densities 
[Eq(25d)]  and  is  a function  of  time  due  to  the  pyrolysis  reactions: 
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The  cell  pyrolysis  for  the  time  step  is  evaluated  at  T = (T„  + T„+j)  / 2,  and  component  densities  also 
at  jS  = 'h.  Eq.(66)  was  rewritten  in  a form  allowing  numerical  solution  by  a standard  secant  method  (ref 
[14],  pp  299-306).  This  completes  the  outline  of  Larkin’s  method.  The  net  result  of  using  Larkin’s 
method  is  that  the  time  step  is  no  longer  limited  by  grid  size,  and  pyrolysis  is  modeled  in  a reasonable 
execution  time. 


4.  EXPERIMENTS 
4a.  Experimental  Arrangement 

The  next  step  is  to  compare  the  results  predicted  by  TMPSUB2  with  experimental  results,  in  which  an 
electrical  heater  was  substituted  for  a cigarette.  The  schematic  of  an  experiment  designed  to  measure  time 
to  ignition  is  shown  in  figure  4.  This  consists  of  the  heater  element  from  an  automobile  cigarette  lighter, 
fitted  with  a concentric  jacket  to  permit  an  air  purge.  The  purpose  of  the  air  purge  is  to  keep  evolved 
products  from  the  sample  from  being  ignited  to  flaming,  rather  than  to  merely  smolder.  The  heating 
element  is  raised  to  varying  temperatures,  all  of  them  high  enough  for  the  element  to  be  glowing  (500  - 
900°C).  The  purging  jet  comes  through  the  jacket,  past  the  face  of  the  element,  and  then  out  normal 
to  it,  toward  the  sample;  it  picks  up  a good  deal  of  heat  as  it  travels  through  the  device.  The  resulting 
flux  distribution  beneath  the  heating  element  is  shown  in  figure  5.  As  we  would  expect,  it  is  axially 
symmetric;  it  is  well  fitted  by  a Gaussian  profile.  Note  that  the  higher  the  temperature  the  disk  is  heated 
to,  the  higher  is  the  peak  flux. 

The  measured  flux  to  the  gauge  is  not  the  same  as  what  the  substrate  sees,  however,  because  the  latter 
heats  up,  whereas  the  flux  gauge  does  not;  call  the  former  0 and  the  latter,  0^.  It  follows  from  Eq(3) 
that 

4>,  = (68) 

where  Tg  is  the  temperature  of  the  hot  purging  gas,  T^.  that  of  the  (cold)  gauge,  and  (f)^  = For  the 
case  where  the  substrate  is  being  heated,  however,  the  convective  contribution  goes  down: 

<!’.  = (69) 

where  T5(t)  is  the  temperature  of  the  substrate,  which  increases  continuously.  Thus  the  total  flux 
decreases  monotonically.  Although  we  do  not  know  the  gas  temperature  Tg,  we  do  not  need  it;  for,  we 
can  combine  Eqs.(68)  and  (69)  to  obtain 

(70) 


The  heat  transfer  coefficient  h was  discussed  in  Section  2.  It  is  found  for  this  experimental  configuration 
in  Appendix  D.  Analysis  of  the  experimental  results,  given  in  Appendix  D,  yields  a number  of  values 
important  for  this  experiment,  including  h and  the  disk  temperatures;  the  results  are  shown  in  Table  Dl; 
h is  found  to  be  a weak  function  of 

When  a mock-up  consisting  of  flexible  PU  foam  covered  by  #12  cotton  duck  was  placed  in  position  under 
the  heat  source,  smoldering  ignition  occurred  after  a certain  amount  of  time;  the  ignition  delay  depends 
on  the  intensity  of  the  flux,  as  indicated  qualitatively  by  Eq.(19),  Section  2b.  The  experimentally- 
obtained  ignition  delays  are  plotted  in  figure  6 as  a function  of  the  peak  heat  flux. 
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Analysis  of  the  ignition  delay  data: 

The  ignition  delay  times  as  a function  of  the  external  flux  are  given  in  Table  1. 


Table  1.  Ignition  delay  times 


«,„(kW/ni2)  - 18 

25 

34.5 

44 

tig(sec)  - 472 

70 

37 

22 

Since  is  constant  for  each  run,  then  according  to  Eq(19),  a plot  of  versus  should  yield  a 
straight  line  whose  intercept  is  the  critical  flux,  When  this  is  done,  however,  the  result  is  not  a 

straight  line.  More  generally,  then,  we  assume  that  the  relationship  has  the  form 


Ad)  ^ ^ 


(71) 


Only  the  correct  choice  for  <})^^■^^  will  yield  a straight  line  in  a logarithmic  plot  of  A<j)  vs  tjg.  When  that 
is  obtained,  the  slope  of  the  line  is  -p;  it  was  found  that  d>grjt  = 16.9  kW/m^  and  that  p = 1.087,  with 
A varying  from  757  to  816;  the  average  value  is  798.  If  we  take  p = 1,  A increases  slowly  with 

The  fact  that  p = 1,  rather  than  1/2,  indicates  that  the  substrate  behaves  like  a thermally  thin  material; 
that,  in  turn,  suggests  that  it  is  principally  the  fabric  which  is  involved  in  the  heating  and  ignition.  If  we 
assume  that  the  principal  cooling  mechanism  is  radiation,  then  <|)^^■^^  = 16.9  kW/m^  implies  that  if  the 
radiative  absorption  coefficient  of  the  fabric  is  a = 0.9,  then  T^g  = 759  K = 485 °C  --  a remarkably 
high  value.  Even  with  a = 1.0,  Tjg  = 739  K = 466°C,  considerably  greater  than  the  400°C  measured 
independently.  This  shows  that  there  is  substantial  convective  cooling  during  the  heating/ignition  process. 


4b.  Material  Data 

In  order  to  see  whether  the  program  can  calculate  tjg  correctly,  it  is  necessary  to  calculate  the  surface 
temperature  under  the  heater.  In  order  to  do  that,  however,  it  is  necessary  to  have  correct  input  data; 
among  these  data  are  the  thermophysical  data,  icpc. 

4b.  1 Fabric 

First,  consider  the  fabric.  The  fabric  that  was  used  in  the  experiment  was  #12  cotton  duck.  Material 
data  for  cotton,  especially  as  a function  of  temperature,  are  surprisingly  difficult  to  find,  even  though  it 
is  a common  and  long-used  material. 

In  Appendix  E,  we  find  that  reasonable  values  for  the  thermophysical  constants  for  this  particular  cotton 
fabric  (#12  cotton  duck)  are 


p = 620  kglm^ 


K(r)  = 0.28505  k/T)  + 0.84554  K^^,(r)  W/m-K  (72) 


25 


where 


kAD  = k,(TJ(TITJ  = 1.242xl0-’r  WIm-K 


(with  T in  Kelvins)  and 


K — 1 7 • 

''■gas  ^ * ”‘air’ 


with 


a^/f 


, bxlO-^^ 
1 + 


(73) 


(74) 


The  constants  are  a = 2.6464x10'^,  b = 245.4,  and  c = 12.  Then  /c(25°C)  = 0.1435  W/m-K  and 
c(T)  = K(J)/pa  = 7819<c(T);  thus,  c = 1122  J/kg-K  at  300  K. 

4b.2  Char 

The  cotton  decomposes  and  pyrolyses  to  char.  According  to  Parker  [15],  [16],  the  specific  heat  of  char 
is  just  about  that  of  carbon: 

c^(T)  = 1.43  + 3.55xl0-^T  - 1 32x10^ IT^  J/g-K,  with  T in  K (75) 

The  thermal  diffusivity  of  wood  char  is  approximately  constant: 

~ 2.1x10'^  m^/s  (76) 

The  density  of  the  char  depends  on  whether  the  fibers  contract  while  pyrolyzing,  or  not.  Finally,  the 
thermal  conductivity  of  char  is  found  from 

/Cc(T)  = a^PcCcCT).  (77) 


4b. 3 Reaction  kinetic  parameters 


When  the  analysis  of  a pyrolyzing  material  is  carried  out,  it  is  done  for  a thin  layer,  in  order  to  ensure 
uniformity  of  temperature  and  of  oxygen  concentration  throughout  the  sample.  It  is  then  simplest  to  find 
the  reaction  rate  in  the  form 


= 


djWJWJ 

dt 


= A, 


exp(-EJRT)  min 


-1 


0/ 


(78) 


for  the  degradation  reaction,  rather  than  as  in  Eq(20);  here  W^j  = Wjj(t)  is  the  (instantaneous)  weight  of 
the  remaining  (virgin)  material  = weight  of  sample  minus  weight  of  char  and  ash,  and  is  the  original 
weight.  Note  that  A is  commonly  given  in  reactions/minute.  The  actual  reaction  rate  is  given  by 

“ Po^d- 

The  relationship  between  (20)  and  (78)  is  simple:  since  and  p^  = Wq/Vq,  then  as  long  as 

Vo/Vd  = const,  WjAVo  = Pd/Po,  and,  as  in  Eq.(20), 
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where 


Pd  exp(-...) 


(in  the  appropriate  units).  T.J.  Ohlemiller  (private  communication,  10/91)  found  that  for  pyrolysis  in 
pure  nitrogen, 

= 7.49x10*'*  min"*,  n = 0.60,  and  = 43,600  cal/mole  = 182.4  kJ /mole 

Instead  of  these  parameters,  however,  we  shall  use  the  results  of  Kashiwagi  and  Nambu  [6]  for  a 
cellulosic  paper  pyrolyzed  in  air;  that  will  maintain  internal  consistency  in  the  reaction  set.  The  global 
kinetic  constants  given  by  Kashiwagi  and  Nambu  for  the  thermal  degradation  reaction  of  their  cellulosic 
paper  are 

Ejj  = 220  kJ/mole  n^,  = 1.8 

Aj  = 1.2x10*^  min"*  = - 570  J/g 

The  parameters  for  the  other  two  reactions  are  given  in  their  Tables  1 and  2,  reproduced  below.  For  the 
oxidative  pyrolysis  reaction. 


= 

op  ox 


W. 


'/op 


exp 


0/ 


op 

RT 


(79) 


where  is  the  volume  fraction  of  oxygen.  The  kinetic  constants  are 

Eop  = 160  kJ/mol  n^jp  = 0.5 

Aop  = 1.5x10*'*  min"*  nfop=  1.3 

and  H^p  = 5.7  kJ/g 

For  the  char  oxidation  reaction. 


with  the  kinetic  constants 
E^  =160  kJ/mol 
Aj.  = 3.4x10**  min"* 
and 

The  fraction  of  ash  which  remains  is  9%  by  mass.  Eqs.(79)  and  (80)  transform  to  the  standard  form  in 
the  same  way  that  (78)  does.  Thus 


( W^] 

»c 

— 

exp 

IF 

\ oj 

[ RT) 

(80) 


nco  = 0.78 
n^  =1.0 
H,  = 25  kJ/g 


^op  ~ ^obP 


1-n, 


'/op  , 


op' 
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4b.4  Heats  of  combustion 


Brandrup  and  Immergut  (ref  [17])  give  the  heat  of  combustion  from  several  measurements;  as  expected, 
the  results  vary: 

Hj. (cotton)  = 18850  J/g 
Hc(fabricl)  « 15450  J/g 

Hc(fabric2)  « 16700  ± 250  J/g  (fabric  weight  was  180  give?) 

Note:  as  material  pyrolyzes,  we  should  take  c(T)  and  x(T)  for  the  combination  of  materials;  for  the  sake 
of  simplicity,  however,  we  assume  that  the  virgin  fuel,  the  pyrolyzed  material,  the  char,  and  even  the 
ash,  all  have  the  same  specific  heat.  Then  the  appropriate  density  to  use  in  Eq.(12),  for  example,  is  the 
density  of  the  solid,  p^,  which  is  the  sum  of  the  bulk  densities  of  all  the  components  of  the  solid  (see 
Eq(25d)).  As  pyrolysis  and  combustion  proceed,  p^  decreases  monotonically  (assuming  no  fabric 
shrinkage). 

4b. 5a  Foam 

For  the  padding  component  of  the  substrate,  we  have  PU  foam,  for  which  we  have  the  following  data 
(T.J.  Ohlemiller,  private  communication,  12/91): 


p = 0.032  g/cm^  = 32  kg/m^ 

(81a) 

Cp  = 1.46J/g-K 

(81b) 

ic(T)  = 0.03613  + 2.003x10-^  T, 

(81c) 

with  T in  °C. 

There  is  one  other  datum  for  PU  foam:  < xpc> , the  mean  value  of  the  thermal  inertia,  was  measured 
in  the  LIFT  apparatus  (see  [5],  [18],  and  [19])  by  Quintiere  and  Harkleroad  [19].  However,  the  foam 
melts  before  it  ignites,  so  that  p (and,  consequently,  k)  go  up  substantially;  moreover,  the  measurement 
entails  all  the  other  model  approximations.  Hence  the  measurement  is  not  directly  useful. 


4b. 5b  PU  Foam  kinetic  parameters 


The  first  pyrolysis  reaction,  which  results  in  the  collapse  of  the  foam  structure,  yields  0.7  g of  liquid  and 
0.3  g of  vapor,  for  every  gram  of  foam  which  pyrolyzes.  The  reaction  rate,  RRj  = (mass 
gasified/min)/(mass  of  foam),  is  given  by 


RR^  = 0.30fi 


[ P - P^ej 

m 

^Po  “ P«y 

exp 

i'/erj 

(82) 


where  m = 1.5  Also, 

B = 3.36x10*^  and 
Eg  = 44,700  cal/mole. 


Pq  = starting  density, 

p = current  density,  and 

Pg  = density  of  residue  = 0.70po 


These  densities  are  those  of  the  solid  phase  (p^  of  order  Ig/cm^),  rather  than  that  of  the  foam  (p^  = 0.03 
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g/cm^).  The  total  heats  of  pyrolysis  of  PU  foams  range  between  100  and  200  cal/g;  for  this  foam,  it  is 
closer  to  the  upper  value.  Assuming  it  varies  linearly  with  the  mass  loss,  we  can  estimate  30%  of  200, 
or 

Hp  ~ 60  cal/g  * 240  J/g 

Finally,  for  the  diffusion  coefficient  for  O2  in  the  fabric,  see  Eq.(28). 

5.  RESULTS  AND  DISCUSSION 

Rather  than  give  the  final  results  immediately,  it  was  decided  to  describe  the  evolution  of  our  thinking, 
including  some  of  the  false  steps  we  took,  because  a good  deal  can  learned  that  way,  and  it  could  be 
instructive  for  anyone  else  who  might  work  on  this  problem  in  the  future.  With  that  in  mind,  we  first 
consider  the  results  of  calculations  made  assuming  an  inert  substrate: 

5a.  Inert  Case 

First,  the  computer  program  was  used  to  calculate  the  peak  surface  temperature  under  the  heater  as  a 
function  of  time,  for  several  peak  fluxes.  The  results  are  shown  in  figure  7. 

The  ignition  temperature  of  cotton  was  measured  to  be  390-400°C  (T.J.  Ohlemiller  et  al,  private 
communication).  According  to  figure  7,  if  Tjg  = 400°C,  then  ignition  is  attained  for  the  44  kW/m* 
exposure  at  t » 18  sec,  and  for  the  34.5  kW/m^  exposure  at  t = 50  sec.  These  values  are  to  be 
compared  with  experimental  values  of  22  and  37  sec,  respectively;  this  is  fairly  good  agreement.  On  the 
other  hand,  ignition  is  not  attained  at  all  for  the  25  and  18  kW/m^  cases;  a conceivable  reason  might  be 
that  the  ignition  temperature  given  above  was  overestimated.  In  order  to  intersect  the  </>  = 18  kW/m^ 
curve  at  t » 472  sec,  the  ignition  temperature  would  have  to  have  been  about  290°C;  it  is  very  unlikely 
that  such  a large  error  in  the  measurement  of  Tjg  would  have  been  made.  Even  if  it  were,  that  would 
yield  ignition  times  of  7,  12,  and  25  sec,  respectively  --  much  shorter  than  the  measured  times  of  22,  37, 
and  70  sec. 

Three  possible  reasons  for  the  calculated  ignition  times  being  so  short  are: 

(a)  Endothermic  pyrolysis  takes  place,  which  slows  the  temperature  rise. 

(b)  There  is  significant  radiative  heat  transfer  within  the  material,  so  that  the  flux  is  absorbed  in  depth, 
rather  than  mainly  at  the  surface. 

(c)  Tjg  = 290°C  and  the  estimates  of  /<pc(T)  were  in  error,  the  actual  value  being  larger. 

Now,  if  we  consider  endothermic  pyrolysis,  then  we  must  also  consider  exothermic  pyrolysis.  Second, 
the  effect  of  radiative  heat  transfer  will  be  small  when  T is  near  T^,  and  small  in  comparison  with  the 
effects  of  pyrolysis  when  T>  >T^.  Finally,  consider  item  (c):  If  /cpc  were  larger,  then  the  rate  of  rise 
of  Tg  would  be  smaller,  as  can  readily  be  seen  from  Eqs(17)-(19).  The  final  temperature,  however,  is 
independent  of  xpc.  Hence  not  only  would  Tjg  have  to  be  much  smaller,  but  the  estimates  of  /cpc  would 
have  to  have  been  9 times  too  small.  This  is  all  possible,  but  extremely  unlikely. 

5b.  Results  with  pyrolysis 

On  the  other  hand,  we  know  that  pyrolysis  must  take  place;  Figure  8 shows  the  effect  of  including 
pyrolysis,  for  the  Q = 25  kW/m^  case.  Curve  A indicates  what  the  temperature  of  the  fabric  surface 
would  be  if  the  fabric  (and  foam)  were  inert.  Curve  B is  the  result  of  "turning  on"  the  (endothermic) 
degradation  reaction.  Note  that  although  the  rate  of  growth  of  temperature  is  slowed  down  quite 
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perceptibly,  as  expected,  the  temperature  the  surface  reaches  asymptotically  is  exactly  the  same  (assuming 
the  same  surface  absorptivity /emissivity  as  before). 

However,  the  exothermicity  of  the  oxidative  pyrolysis  and  of  char  oxidation  are  far  greater  than  the 
endothermicity  of  thermal  decomposition,  and  should  overwhelm  it.  We  cannot,  therefore,  consider 
thermal  decomposition  alone:  if  we  consider  pyrolysis,  we  must  include  ^ the  major  reactions.  Another 
way  of  seeing  this  is  to  note  that  in  order  to  have  ignition,  we  must  have  a "thermal  runaway,"  where 
exothermic  pyrolysis  heats  up  the  material  faster  than  heat  diffusion  and  surface  losses  can  carry  the  heat 
away. 

Curve  c indicates  what  takes  place  when  the  oxidative  pyrolysis  reaction  is  included,  as  well:  it  is 
exothermic,  and  adds  a great  deal  more  energy  than  the  degradation  reaction  removes.  The  result  is  that 
the  temperature  rises,  rather  than  falling.  Indeed,  it  rises  so  rapidly  that  it  begins  to  appear  like  a thermal 
run-away.  However,  the  temperature  reaches  a peak,  than  declines.  The  reason  is  simple:  the  fuel  is 
rapidly  exhausted;  once  that  happens,  the  heat  source  is  reduced  to  the  original  external  flux.  This  is 
clearly  depicted  in  Figure  9,  which  shows  the  fuel  density  as  a function  of  time:  The  ordinate  is 
temperature,  in  °C,  md  density  in  kg/m^.  The  curves  marked  Tq  and  pq  correspond  to  the  cell  with  the 
highest  temperature.  Those  marked  Tj  and  pj  correspond  to  the  (laterally)  adjacent  cells,  and  those 
marked  T2  and  P2  to  the  next  ring  of  cells.  These  calculations  were  carried  out  assuming  no  char 
oxidation.  We  see  that  as  the  reaction  accelerates,  T^  "runs  away"  and  the  density  plummets  from  its 
virgin  value  to  that  of  char.  Most  significantly,  the  peak  temperatures  develop  just  after  the  density  falls. 

Finally,  curve  d of  Figure  8 shows  the  result  of  adding  the  char-oxidation  reaction  as  well:  that  additional 
heat  source  takes  the  pyrolysis  "over  the  top":  the  temperature  continues  to  run  away  (we  have  arbitrarily 
cut  off  the  calculation  at  600°C,  here).  The  question  then  arose:  why  is  (was)  there  an  oscillation  in  T^? 
The  explanation  is  qualitatively  clear:  if  a cell  is  (too)  large,  then  the  surface/volume  ratio  is  small,  and 
the  heat  cannot  diffuse  away  rapidly  enough.  This  is  confirmed  by  Figure  10:  with  a cell  size  taken  to 
be  a bit  less  than  half  as  large,  the  amplitude  of  the  oscillation  declined  considerably,  and  with  the  cell 
size  halved  again,  the  oscillations  have  almost  disappeared. 

Presumably,  as  we  continue  to  decrease  the  cell  size,  the  numerical  errors  should  become  progressively 
smaller,  until  a further  decrease  in  cell  size  would  make  no  further  difference  in  the  results.  However, 
it  was  found  that  the  results  were  apparently  not  converging  with  decreasing  Ax,  and  that  the  sensitivity 
varied  with  Q.  Investigation  of  the  reason  for  this  great  sensitivity  revealed  that  it  lay  in  our  use  of  a 
variable  grid  size:  whereas  for  grids  of  constant  spacing,  the  numerical  approximations  are  correct  to 
second  order  in  Ax,  that  accuracy  drops  to  somewhere  between  first  and  second  order.  Indeed,  if  the 
numerical  error  is  proportional  to  (Ax)",  then  n is  given  by  a complex  formula  (see  ref  [20],  pp.43  and 
51).  The  "flavor"  of  that  expression  is  given  by  n ~ (2-l-a)/(l-l-a),  where  a = dAxIdx.  That  is,  the 
accuracy  depends  inversely  on  the  rate  of  growth  of  the  grid  size.  An  explicit  calculation  is  given  in 
Appendix  A;  see  Table  A1  there. 

By  reducing  the  rate  of  increase  of  grid  size,  therefore,  the  accuracy  was  increased,  and  the  results  made 
to  converge  better.  Moreover,  it  was  decided  to  switch  from  an  explicit  solution  method  to  Larkin’s 
semi-implicit  method  (see  Section  3).  After  those  changes,  we  arrived  at  the  results  shown: 

The  results  of  four  calculations  are  shown  in  Figure  1 1 : these  were  made  with  an  assumed  heat  transfer 
coefficient  h = 20,  and  an  assumed  oxygen  concentration  of  20%.  For  initial  fluxes  with  peak  values 
of  25  and  34.5  kW/m^,  calculations  are  first  made  assuming  that  the  substrate  is  inert,  and  then  pyrolysis 
is  "turned  on."  We  see  that  for  the  34  kW  case,  the  asymptotic  temperature  lies  at  about  430°C,  and  that 
a thermal  runaway  begins  at  about  300°C,  at  t » 20,  and  is  completed  at  t = 25  sec;  that  would  then 
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be  the  ignition  time.  It  is  difficult  to  state  precisely  what  "the  ignition  temperature"  is,  in  this  case  (but 
see  below). 

For  the  25  kW  case,  the  asymptotic  temperature  lies  at  about  360 °C;  the  first  perceptible  deviation  of 
the  curve  for  the  reactive  case  from  that  for  the  inert  material  occurs  at  T = 280°C,  and  is  clearly 
established  at  300°C.  The  thermal  runaway  takes  place  at  t = 48  sec.  These  times  may  be  compared 
with  the  experimentally-established  ignition  times,  which  were  about  22  and  70  sec,  respectively.  Thus 
we  have  achieved  semi-quantitative  agreement.  In  Figure  12,  on  the  other  hand,  we  see  that  a thermal 
runaway  --  and  therefore  ignition  --  apparently  did  not  take  place  for  the  18  kW  case.  Here,  the 
asymptotic  temperature  is  about  290 °C;  this  is  apparently  not  high  enough  to  produce  a runaway. 

Wherein  lies  the  difficulty?  Most  likely,  one  or  more  of  the  input  values  is  incorrect.  Altering  xpc 
would  principally  change  the  time  scale.  Fig.  13  shows  the  effect  of  changing  the  assumed  oxygen 
concentration:  the  highest  curve  reproduces  the  upper  curve  in  Fig.  12.  The  next  two  curves  show  what 
happens  when  < Y>  is  assumed  to  be  0.15  and  0.11,  respectively.  Finally,  the  bottom  curve  (again) 
corresponds  to  the  inert  case.  Thus,  it  is  not  [O2]  having  been  chosen  too  low  that  prevents  ignition,  and 
one  or  more  of  the  kinetic  parameters  is  probably  in  error.  Refer  to  Figure  14:  curves  a and  b again 
reproduce  Fig.  12  (on  a different  scale).  For  curve  c,  the  char-oxidation  rate  was  doubled.  It  is  apparent 
that  the  curves  overlap  completely.  A preliminary  conclusion  inferred  from  this  was  that  the  observed 
result  was  due  to  all  the  char  that  is  produced  already  being  oxidized.  Observation  of  the  kinetic 
constants  for  oxidative  pyrolysis  and  char  oxidation,  however,  makes  it  clear  that  the  latter  is  three  orders 
of  magnitude  slower  than  the  former.  Therefore  merely  doubling  the  char-oxidation  rate  will  only  perturb 
the  energy  output  slightly  - so  slightly  that  it  will  not  even  show  up  in  the  figure. 

For  curve  d,  the  oxidative  pyrolysis  rate  was  doubled  (the  preexponential  factor  A was  doubled), 
doubling  the  char-production  rate  through  this  branch;  this  indeed  produced  a thermal  runaway.  Curves 
e and  f are  the  results  of  increasing  A by  20%  and  10%,  respectively.  Thus,  a quite  modest  increase  in 
A produces  (predicts,  that  is)  ignition,  although  at  230  sec,  rather  that  the  measured  472  sec).  Such  an 
increment  is  not  only  well  within  experimental  error,  but  --  more  to  the  point  - is  entirely  plausible, 
when  the  likely  differences  between  the  cellulosic  paper  and  a cotton  fabric  (with  different  impurities)  are 
considered.  On  the  other  hand,  it  was  assumed  that  y = 0.20.  For  y = 0.15,  A must  be  increased  by 
30%  in  order  to  get  ignition  (the  resulting  T(t)  curve  is  very  similar  to  the  "best"  one).  Although  this 
is  greater  than  the  10%  increase  found  above,  it  is  still  entirely  plausible. 

We  have  so  far  considered  the  sensitivity  of  the  results  to  the  oxygen  concentration,  the  thermophysical 
constants  of  the  fabric,  and  its  kinetic  parameters.  Surprisingly,  there  are  two  other  significant 
parameters:  first,  if  we  use  h = 22  for  the  heat  transfer  coefficient  in  Eq.(69)  (as  suggested  by  the  results 
of  Appendix  D,  shown  in  Table  Dl),  rather  than  the  assumed  h = 20,  the  "asymptotic"  temperature  (that 
at  t = 500  sec)  reaches  only  272 °C,  rather  than  292  °C,  and  there  is  not  the  faintest  possibility  of 
achieving  ignition,  unless  is  substantially  smaller  than  160  kJ/gm.  One  solution  is  to  assume  that 
h=20  is  the  correct  value  to  use,  since  the  theoretical  calculations  in  Appendix  D could  easily  be  off  by 
10%  or  more.  Another  resolution  is  possible,  too:  the  second  parameter  which  is  important  in  this 
threshold  region  ("threshold,"  because  18  kW/m^  is  close  to  the  critical  flux,  16.9  kW/m^)  is  the  thermal 
conductivity  of  the  foam  padding.  In  all  the  runs  made  above,  it  was  assumed  that  k = 0.056  J/m-K  and 
c = 1.9  J/g-K,  for  the  foam  at  T = 20°C.  If  it  were  assumed  that  k = 0.096,  instead,  then  the  foam 
would  act  as  a more  efficient  heat  sink,  and  the  surface  temperature  could  be  expected  to  drop;  indeed, 
a calculation  showed  that  the  peak  surface  temperature  at  t = 500  s fell  to  272 °C  for  the  inert  fabric. 
In  fact,  however,  the  value  0.056  for  the  thermal  conductivity  was  for  a foam  of  density  48  kg/m^! 
Transforming  that  for  a 32  kg/m^  foam,  according  to  Eq.(62),  yields  0.036,  almost  exactly  what  is  given 
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in  the  present  test  (see  Eq.(81c)).  This  reduces  the  heat  sink,  and  must  yield  an  increased  asymptotic 
surface  temperature.  Using  the  foam  parameters  given  by  Eqs.(81),  and  h = 20,  the  t=500  temperature 
indeed  rises,  from  292 °C  to  300°C;  with  h = 22,  it  is  294°C,  and  we  only  need  to  increase  A by  10% 
to  get  runaway. 

The  "best"  set  of  parameters,  then,  is  that  given  in  sections  4b.  1 and  4b. 3 for  the  fabric,  and  4b. 5a  for 
the  foam.  For  the  heat  transfer  coefficient,  use  the  values  in  Table  Dl.  With  that  set,  we  obtain  the 
curves  shown  in  Fig.  15.  The  corresponding  calculated  ignition  times  are  given  in  Table  2.  Thus,  the 
calculated  values  are  all  about  half  the  measured  values. 

Note  that  the  polyurethane  foam  begins  to  melt  and  recede  from  the  fabric  when  its  temperature  reaches 
about  300 °C,  thereby  decreasing  the  heat-sink  effect  of  the  padding,  and  accelerating  the  heating  of  the 
fabric;  this  effect  has  not  been  included  in  the  model,  however. 


Table  2.  Calculated  vs  experimental  ignition  delay  times  for  the  four  fluxes 

Q (kW/m^)  Ignition  delay,  tjg  (in  sec) 

Calculated  Measured 


18 

209 

472 

25 

37.2 

70 

34.5 

19.0 

37 

44 

12.6 

22 

It  has  been  suggested  that  we  might  avoid  the  necessity  of  explicitly  including  the  pyrolysis  reactions  by 
choosing  some  effective  ignition  temperature.  This  is  not,  in  fact,  feasible:  if  we  take  the  measured 
ignition  times  and  mark  them  on  the  four  curves  in  Figure  7,  corresponding  to  the  inert  assumption,  we 
find  that  they  intersect  these  curves  at  widely  varying  temperatures;  see  Table  3.  It  is  apparent  that  this 
"ignition  temperature"  is  a strong  function  of  the  external  flux. 


Table  3.  Surface  temperatures  which  would  be 
attained  by  the  substrate  at  the  measured  ignition 
times  if  the  substrate  had  been  inert 


(in  kW/m^)  = 

18 

25 

34.5 

44 

T 

•g 

(in  °C) 

291 

342 

380 

400 

T 

•g 

(in  “K) 

564 

615 

653 

673 

On  the  other  hand,  if  we  define  "the  ignition  temperature"  as  the  point  on  the  Tg(t)  curve  where  the 
temperature  is  rising  at  some  rapid  rate  - say,  100°C/sec  - then  from  Figs.  11  and  13  we  see  that  that 
gradient  is  attained  at  the  approximate  temperatures  shown  here: 
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(in  kW/m^)  = 

18 

25 

34.5 

Tn,ax  (in  °C)  = 

390 

390 

380 

T.ax  (in  °K)  = 

663 

663 

653 

These  calculated  temperatures  are 

close 

to  each  other,  and  close  to  the  390-400°C  which  has  been 

measured. 


6.  SUMMARY 

We  have  created  several  computer  programs,  of  which  the  central  one  is  TMPSUB2.  This  program 
calculates  the  temperature  history  at  every  point  in  a substrate  which  is  subjected  to  a strongly  localized 
heating  flux  on  its  top  surface.  The  (solid)  substrate  consists  of  two  layers,  the  top  one  being  a fabric 
and  the  lower  a foam  pad;  there  may  be  a thin  intervening  air  gap.  The  substrate  is  taken  to  be  a 
rectangular  parallelepiped,  and  it  is  broken  up  into  several  thousand  cells.  There  is  a user-friendly  front 
end  for  the  input,  described  in  Section  3.  The  program  runs  well  on  a 386-level  computer  with  a math 
coprocessor,  or  a 486-chip  computer. 

As  stated  in  the  introduction,  this  program  serves  to  calculate  the  temperature  of  the  upholstered  furniture 
as  a function  of  time  and  position,  when  it  is  exposed  to  a prescribed  heating  flux.  This  flux  can  be 
highly  peaked  at  a point,  vary  with  time,  and  move  at  a constant  (specified)  rate  over  the  top  surface  of 
the  furniture,  assumed  to  be  horizontal.  The  radiative  and  convective  heat  losses  from  the  surface  are 
given  correctly.  If  and  when  the  temperature  rate  of  rise  at  a given  location  suddenly  "accelerates"  to 
a value  high  enough  that  the  surface  glows  (that  is,  T > 5()0°C  or  so),  we  can  say  that  smoldering 
ignition  has  occurred.  The  ambient  oxygen  level  can  be  set  at  whatever  value  one  wishes.  The  program 
will  not  tell  us  whether  flaming  ignition  takes  place.  It  also  does  not  treat  the  case  where  the  flux  is 
applied  in  a crevice,  such  as  is  formed  between  the  seat  cushion  and  the  seat  back.  The  program  does 
not  take  into  account  the  effect(s)  of  the  foam  possibly  melting  and  receding  away  from  the  fabric. 
Finally,  it  also  does  not  take  oxygen  diffusion  within  the  cushion  explicitly  into  account;  hence  in  certain 
threshold  situations,  where  a small  change  in  oxygen  concentration  determines  whether  ignition  does  or 
does  not  take  place,  the  results  are  ambiguous  and  not  to  be  trusted.  Note  that  it  is  often  difficult  to 
obtain  the  needed  kinetic  and/or  thermophysical  parameters  for  the  material;  or,  when  available,  to  know 
how  accurately  they  are  known.  Therefore  this  caveat  must  also  be  made:  even  if  the  program  were 
perfect,  its  results  are  only  as  good  as  the  input  parameters  which  are  supplied.  On  the  other  hand,  it 
should  accurately  reproduce  (or  predict)  trends. 

When  the  flux  to  which  the  substrate  is  exposed  is  near  the  critical  flux,  the  result  (i.e.,  ignition  does  or 
does  not  take  place)  is  sensitive  to  most  of  the  parameters,  such  as  the  ambient  oxygen  density,  the  heat 
transfer  coefficient,  the  thermal  conductivity  of  the  foam  as  well  as  that  of  the  fabric,  etc. 

The  way  T(r,t)  is  found  is  by  solving  the  PDE  which  describes  the  diffusion  of  heat  in  a solid,  Eq.(12), 
numerically.  The  solid  is  subjected  to  a nonuniform  and  time-varying  heating  flux  at  its  top  surface,  and 
(simultaneously)  experiences  convective  and  radiative  heat  losses.  Moreover,  the  solid  can  undergo 
pyrolitic  reactions;  we  consider  three,  here:  one  endothermic  step  (thermal  degradation  to  char),  one 
oxidative  pyrolysis  to  char,  and  oxidation  of  the  char  (to  ash). 
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The  equation  set  is  very  stiff,  because  of  the  highly  nonlinear  form  of  the  (Arrhenius)  reactions.  We  have 
therefore  used  a semi-implicit  method  to  solve  the  equation  set. 

The  way  ignition  is  seen  to  be  achieved  (in  a cell)  is  that  the  temperature  undergoes  a "thermal  runaway." 
This  does  not  occur  abruptly  as  a particular  temperature  is  reached,  so  that  the  "ignition  temperature"  is 
not  well  defined.  When  the  transition  region  between  "inert"  heating  and  thermal  runaway  is  narrow, 
the  concept  is  adequate.  According  to  this  model,  the  transition  region  is  rather  broad,  and  so  it  is  not 
so  useful  a concept  in  this  context. 

What  is  and  is  not  in  the  program  is  listed  in  Table  1,  section  2. 

Note  that 

a)  We  can  get  a first  approximation  to  the  radiative  heat  transfer  within  the  material  by  using  Kunii’s 
expression,  Eq.(14). 

b)  The  effects  of  cation  concentration  can  probably  be  modeled  by  appropriate  changes  (not  described 
here)  in  the  kinetic  parameters.  Finally, 

c)  the  effect(s)  of  relative  humidity  in  the  ambient  can  likewise  probably  also  be  approximately  modeled 
by  making  appropriate  changes  (again,  not  described  here)  in  p,  c,  and  k. 

An  experiment  was  carried  out  to  ignite  the  substrate.  In  trying  to  reproduce  those  experimental  results, 
it  was  found  that  the  calculated  results  are  sensitive  to  the  input  values  chosen,  especially  the  kinetic 
parameters.  It  was  found  that  the  preexponential  factor  found  by  Kashiwagi  and  Nambu  (ref  [6])  for  the 
oxidative  pyrolysis  reaction  in  a cellulosic  paper  had  to  be  increased  by  10%  for  the  cotton  duck  fabric 
in  order  to  get  ignition  for  the  lowest  of  a set  of  heating  fluxes  to  which  the  substrate  was  exposed.  The 
result  was  semi-quantitative  agreement  with  the  observed  ignition  times. 
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FIGURE  CAPTIONS 


Figure  1.  One-dimensional  conduction 

Figure  2.  Air  gap  heat  transfer  model 

Figure  3.  Substrate  coordinate  system 

Figure  4.  Schematic  of  heat  source  for  ignition  tests 

Figure  5.  Flux  profile  from  heat  source,  as  measured  by  a total  heat  flux  gauge 

Figure  6.  Time  required  to  ignite  the  fabric,  for  different  (initial)  heat  flux  exposures.  Crosses 
correspond  to  the  discussion  in  the  text.  Filled  circles  correspond  to  a different  set  of 
experiments. 

Figure  7.  Peak  surface  temperatures  of  substrate  as  a function  of  time,  for  the  four  exposures 

Figure  8.  Peak  surface  temperatures  of  substrate  exposed  to  Q = 25  kW/m^,  as  a function  of  t,  for 

a.  Material  assumed  to  be  inert  d.  Char  oxidation  also  included,  as  well 

b.  Thermal  degradation  only  e.  Only  one  reaction:  oxidative  pyrolysis 

c.  Oxidative  pyrolysis  as  well 

Figure  9.  Temperature  and  density  of  central  cell  (subscript  0),  adjacent  cells  (subscript  1),  and  cells  in 
next  ring  around  center  (subscript  2),  as  functions  of  time 

Figure  10.  Temperature  of  central  surface  cell  (i.e.,  peak  temperature)  for  three  different  grid  sizes:  1 .25, 
0.50,  and  0.25-mm  cubes.  Char  oxidation  was  purposely  left  out 

Figure  11.  Peak  temperature  for  the  25  and  34.5  kW/m^  cases,  assuming  (a)  no  pyrolysis,  and  (b)  all 
three  pyrolytic  reactions 

Figure  12.  Peak  temperature  for  the  18  kW/m^  case,  with  the  same  assumptions  as  in  Fig.  11 

Figure  13.  Peak  temperature  for  the  18  kW/m^  case,  for  several  values  of  mean  O2  mass  fraction,  < y > 

Figure  14.  Peak  temperature  for  the  18  kW/m^  case,  with  various  assumptions  for  the  pyrolysis: 

Curve  a,  no  pyrolysis;  curve  b,  "standard"  pyrolysis;  curve  c,  double  the  char  oxidation  rate; 
Curve  d,  double  the  oxidative  pyrolysis  rate;  Curve  e,  1.2  times  the  oxidative  pyrolysis  rate; 
Curve  f,  1 . 1 times  the  oxidative  pyrolysis  rate 

Figure  15.  Peak  temperature  as  a function  of  time,  for  all  four  cases,  using  the  best  set  of  input  data 
Figure  16.  Total  heat  flux  impinging  on  gauge,  with  and  without  purge  flow 

Figure  17.  Calculated  heat  transfer  coefficient,  as  a function  of  the  temperature  Tg  of  the  impinging  purge 
gas  Jet 

Figure  18.  Thermal  conductivity  of  cotton  as  a function  of  temperature,  as  measured  by  different  workers 
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Figure  4.  Schematic  of  heat  source  for  ignition  tests 
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Figure  5.  Flux  profile  from  heat  source,  as  measured  by  a total  heat  flux  gauge 
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Figure  6.  Time  required  to  ignite  the  fabric,  for  different  (initial)  heat  flux  exposures.  Crosses 
correspond  to  the  discussion  in  the  text.  Filled  circles  correspond  to  a different  set  of 
experiments. 
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Figure  7.  Peak  surface  temperatures  of  substrate  as  a funaion  of  time,  for  the  four  exposures 
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Figure  8,  Peak  surface  temperatures  of  substrate  exposed  to  Q = 25  kW/m^,  as  a function  of  t,  for 


a.  Material  assumed  to  be  inert 

b.  Thermal  degradation  only 

c.  Oxidative  pyrolysis  as  well 


d.  Char  oxidation  also  included,  as  well 

e.  Only  ^ reaction:  oxidative  pyrolysis 
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Figure  9.  Temperature  and  density  of  central  cell  (subscript  0),  adjacent  cells  (subscript  1),  and  cells  in 
next  ring  around  center  (subscript  2),  as  functions  of  time 
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Figure  10.  Temperature  of  central  surface  cell  (i.e.,  peak  temperature)  for  three  different  grid  sizes.  1.25, 
0.50,  and  0.25-mm  cubes.  Char  oxidation  was  purposely  left  out 
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Figure  11.  Peak  temperature  for  the  25  and  34.5  kW/m^  cases,  assuming  (a)  no  pyrolysis,  and  (b)  all 
three  pyrolytic  reactions 
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Figure  12.  Peak  temperature  for  the  18  kW/m^  case,  with  the  same  assumptions  as  in  Fig.  11 
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Figure  13.  Peak  temperature  for  the  18  kW/m^  case,  for  several  values  of  mean  O2  mass  fraction,  <y> 
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Figure  14.  Peak  temperature  for  the  18  kW/m^  case,  with  various  assumptions  for  the  pyrolysis; 

Curve  a,  no  pyrolysis;  curve  b,  "standard"  pyrolysis;  curve  c,  double  the  char  oxidation  rate; 
Curve  d,  double  the  oxidative  pyrolysis  rate;  Curve  e,  1.2  times  the  oxidative  pyrolysis  rate; 
Curve  f,  1 . 1 times  the  oxidative  pyrolysis  rate 
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Figure  15.  Peak  temperature  as  a function  of  time,  for  all  four  cases,  using  the  best  set  of  input  data 


49 


EFFECT  OF  PURGE  FLOW  ON  HEAT  FLUX  VS  HEATER  PO\fER 
MEASURED  AT  PEAK  POSITION;  5.4  MM  BELOW  HEAT  SOURCE 
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Figure  16.  Total  heat  flux  impinging  on  gauge,  with  and  without  purge  flow 
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Figure  17.  Calculated  heat  transfer  coefficient,  as  a function  of  the  temperature  Tg  of  the  impinging  purge 
gas  jet 
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Figure  18.  Thermal  conductivity  of  cotton  as  a function  of  temperature,  as  measured  by  different  workers 
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TEMPERATURE 


APPENDIX  A 


Conduction  Algorithm  Tests 


A computer  program  called  TEMPSUB  was  developed  as  part  of  the  earlier  investigation  into  the  ignition 
of  furnishings  by  smoldering  cigarettes  (Gann  et  al,  ref  [2]).  This  program  modeled  heat  transfer  in 
furniture,  or  in  a "substrate,"  using  a simple  finite  difference  approximation  (FDA)  for  a homogeneous 
substrate  with  uniform  and  constant  properties.  The  research  indicated  that  this  program  would  have  to 
be  expanded  to  include  a two-layer  model  (fabric  -I-  padding),  pyrolysis  of  each  layer,  an  asymmetric  flux 
input,  and  a variable  grid. 

These  features  have  been  implemented  by  using  a slightly  different  approach  to  the  FDA  than  was  used 
in  TEMPSUB.  The  original  approach  was  to  convert  the  differential  equation  for  heat  transfer  into  an 
FDA  by  a Taylor’s  series  approximation.  The  new  approach  is  based  on  the  conservation  of  energy 
within  a control  volume.  It  is  based  on  physical  reasoning  and  is  usually  easy  to  apply.  It  is  most  useful 
for  variable  grids,  convective  boundary  conditions,  odd-shaped  regions,  etc.  It  is  more  difficult  to  obtain 
accuracy  estimates  for  the  control  volume  approach  than  for  the  Taylor’s  series  approach.  For  simple 
cases  involving  uniform  grids  and  homogeneous  materials  the  two  approaches  lead  to  identical  FDAs. 
See  Torrance  or  Croft  & Lilley  (refs  [21]  and  [22])  for  further  details. 

Since  there  is  a considerable  increase  in  the  desired  capabilities  of  the  upgraded  program,  it  was  decided 
to  develop  a program  which  would  allow  extensive  testing  of  the  FDA.  This  program  is  called  CTEST3 
(Conduction  TEST  - 3 dimensional)  which  has  the  ability  to  model  simple  boundary  conditions  (constant 
temperature,  heat  flux,  or  convection  coefficient)  on  any  of  the  six  faces  of  the  region. 

The  FDA  used  in  CTEST3,  the  explicit  Euler  method,  has  been  checked  against  several  heat  transfer 
problems  which  have  analytic  solutions.  The  first  few  tests  involve  various  combinations  of  constant 
temperature,  heat  flux,  and  convection  coefficient  boundary  conditions  with  analytic  solutions  from  the 
classic  text  by  Carslaw  and  Jaeger  (ref  [3]). 

Several  of  these  analytic  solutions  involve  the  error  function  which  is  defined  by 


X 


(Al) 


SO  that 


erf(0)  = 0, 


erf(oo)  = 1, 


and 


erf(-x)  = - erf(x). 


The  complementary  error  function  is  also  used.  It  is  defined  as 


(A2) 


so  that 


erfc(0)  = 1, 


and 


erfc(oo)  = 0. 
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Repeated  integrals  of  the  error  function  are  also  useful  in  conduction  problems.  These  are  defined  by 
the  recursive  relationships 


i“erfc(jc)  = J *erfc(jc)<ic  n=l,2,... 

X 

(A3) 

i®erfc(x)  = erfc(x) 

(A3a) 

ierfc  (x)  = i^erfc(x)  = - xerfc(x) 

(A3b) 

2ni“erfc(x)  = i“'^erfc(x)  - 2xi“''erfc(x)  n=2,3,... 

(A3c) 

See  Appendix  II  of  Carslaw  and  Jaeger  [3]  for  further  details  on  error  functions.  Computer  subroutines 
were  written  implementing  these  functions  for  the  computation  of  the  analytic  solutions  of  the  heat 
transfer  tests. 


Test  1:  One-dimensional  steady-state  conduction. 

CTEST3  has  been  tested  for  steady  state  conduction  (with  constant  thermal  properties  and  uniform  grid 
spacing).  This  test  consists  setting  opposite  faces  on  a cubic  region  to  different  temperatures  and  making 
the  remaining  faces  adiabatic.  After  a sufficient  number  of  time  steps  the  temperature  within  the  region 
should  vary  linearly  from  the  hot  to  the  cold  face. 

T\x)  = TXO)  + - [JID  - 7X0)]  (A4) 

L 

This  has  been  confirmed  in  all  three  directions. 


Test  2:  One-dimensional  transient  conduction,  constant  heat  flux  boundary  condition. 


This  test  was  used  in  the  development  of  the  original  substrate  model.  The  analytic  case  for  a constant 
flux  involves  a homogeneous  solid  occupying  the  semi-infinite  region  x > 0.  The  solid  is  initially  at  zero 
temperature  throughout.  At  time  t = 0 a constant  heat  flux,  q,  is  applied  to  the  x = 0 surface.  The 
temperature  within  the  region  is  given  by  (ref  [3],  p 75,  eq  6) 


nx.t)  = 


ierfc 


(A5) 


Preliminary  testing  (again  using  a uniform  grid  and  constant  thermal  properties)  indicated  that  as  At  and 
Ax  decreased,  there  was  a uniform  approach  to  the  analytic  solution.  As  for  the  accuracy  to  be  expected, 
for  Ax  = 1 mm  and  At  = 0.5  sec,  the  error  (after  the  first  three  time  steps)  was  < 1%.  We  note  that 
the  results  of  test  did  not  agree  with  results  from  the  original  TEMPSUB  model.  Further  investigation 
indicated  an  error  in  the  TEMPSUB  boundary  conditions  subroutine.  Correcting  this  error  brought 
results  from  the  two  programs  into  complete  agreement. 
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Tests  show  that  reducing  the  grid  size  along  with  corresponding  reduction  of  the  time  step  cause  the  FDA 
solution  to  approach  the  analytic  solution.  Therefore,  the  FDA  is  consistent.  Reducing  the  time  step 
without  changing  the  grid  size  does  not  improve  the  accuracy  of  the  solution.  In  fact,  it  is  best  to  operate 
as  close  to  the  stability  limit  as  possible  for  both  accuracy  and  execution  time.  Use  of  the  variable  grid 
gives  results  consistent  with  the  uniform  grid  at  the  surface.  The  error  in  the  calculated  surface 
temperature  goes  down  as  time  increases. 


Test  3:  One-dimensional  transient  conduction,  constant  convection  coefficient  boundary  condition. 


This  test  represents  a slab  of  a homogeneous  solid  of  thickness  2L  in  the  x direction  and  infinite  extent 
in  the  y and  z directions  which  is  initially  at  unit  temperature  throughout.  At  time  t=0  the  temperature 
of  the  fluid  on  both  sides  of  the  slab  is  changed  to  zero  and  heat  is  convected  from  the  slab  through  a 
constant  convection  coefficient.  Because  of  symmetry  this  problem  is  equivalent  to  a slab  of  thickness 
L with  one  adiabatic  surface  at  x = 0 and  a convective  surface  at  x = L.  The  temperature  within  the 
region  is  given  by  (ref.[3],  p 122,  Eq(12)) 


2Bcos(6„x/I) 


B=i  [B  +B  + 6„  ]cos6„ 

*■  /I  *■  II 


where  B = hL//c  (Biot  number), 

F = aiH?  (Fourier  number),  and 

are  the  solutions  of  the  transcendental  equation  S^tanS^  = B. 


(A6) 


In  order  to  check  the  calculation  of  this  complicated  analytic  solution,  the  solution  to  a related  problem 
was  also  computed.  This  is  the  temperatures  in  a semi-infinite  slab  with  the  convective  boundary 
condition  (ref. [3],  p 72,  eq(5)) 


II 

n 

( X ) 

- exp 

— (x+a^f/K))erfc 

^ h rz 

— ^ +-y/at 

[yJAat] 

U j 

(A7) 


where  x is  now  the  distance  from  the  convective  surface  into  the  region.  There  is  good  agreement 
between  the  FDA  and  analytic  solutions.  Again,  the  error  in  the  calculated  temperature  goes  down  as 
time  increases. 


Test  4:  Three-dimensional  transient  conduction,  constant  surface  temperature  boundary  condition. 


This  test  represents  a block  of  a homogeneous  solid  in  the  region  defined  by  -a  < x < a,  -b  < y < b, 
and  -c  < z < c which  is  initially  at  unit  temperature  throughout.  At  time  t=0  the  temperatures  of  the 
surfaces  of  the  block  are  reduced  to  zero,  and  the  block  begins  to  cool.  The  temperature  within  the 
region  is  given  by  (ref[3],  p 184,  Eq(5)) 
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T(.x,y,z,t)  = — E E E 


(-1) 


l*m*n 


(2M)(2m+l)(2/i+l) 


cos 


(2/+1)7cx 


2a 


cos 


(2m+l)7ty 


2b 


cos 


(2n+l)nz 


2c 


(A8) 


f 

(2/+l)2  (2w+1)2  (2n+l)2' 

1 

J 
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This  expression  requires  the  summation  of  many  terms  at  small  values  of  time,  but  only  a few  terms  at 
large  values  of  time.  The  implementation  of  this  complicated  equation  had  to  be  checked  against  simpler 
analytic  cases.  The  first  case  represents  a homogeneous  solid  occupying  the  semi-infinite  region,  x > 
0.  The  solid  is  initially  at  unit  temperature  throughout.  At  time  t = 0 the  temperature  at  x = 0 is 
instantly  reduced  to  zero.  The  temperature  within  the  region  is  given  by  (ref[3],  p 59,  eq(3)) 


r(x,f)  = erf 


(A9) 


The  second  case  involves  a solid  which  occupies  the  region  x > 0,  y > 0,  z > 0.  It  is  initially  at  unit 
temperature  and  at  time  t = 0 the  temperature  at  the  x = 0,  y = 0,  and  z = 0 surfaces  is  instantly 
reduced  to  zero.  The  temperature  within  the  region  is  given  by  (refI3],  p 184,  Eq(l)) 


T{x,y,z,t)  = erf 


erf 

( X ] 

erf 

y ] 

erf 

r z 

[p^tj 

[v/4afj 

(AlO) 


The  temperatures  near  the  comers  of  the  block  should  be  very  similar  to  this. 

There  was  good  agreement  (error  < 1 %)  for  a test  with  a=30  mm,  b=20  mm,  c=  10  mm  using  a 1 mm 
uniform  grid.  The  original  variable  grid  model  was  found  to  be  insufficiently  accurate  at  points  where 
the  grid  size  changed.  It  was  therefore  replaced  by  the  current  uniformly  increasing  grid  at  a cost  of 
some  increase  in  code  complexity;  although  the  results  are  not  quite  as  accurate  as  for  the  uniformly- 
spaced  grid,  the  difference  is  very  minor  (see  Section  3d). 


Test  5:  One-dimensional  transient  conduction,  two  different  materials. 


This  test  consists  of  one  material  in  the  region  1 (0  < z < L)  initially  at  unit  temperature  and  another 
material  in  region  2 (z  > L)  initially  at  zero  temperature.  The  boundary  at  z = 0 is  adiabatic.  At  time 
t=0  heat  begins  to  be  conducted  between  the  two  regions.  The  temperatures  in  the  two  regions  are  given 
by  (Ozisik,  ref[23],  p.328,  Eq  8-109) 


where 


jd"  Jerfc 

2 iis  ' 


{ln+\)L-x 


+ erfc 


(2n+l)L+x 


2 B-o 


2nL+\L{x-L) 


- erfc 


(2n+2)L+\x(x-t) 


(All) 


1*  = 


— , P = 


6=1^ 

P^l 


There  was  good  agreement  between  the  FDA  and  analytic  solutions.  It  even  worked  well  when  the  first 
layer  was  only  one-half  of  a grid  thick.  This  may  be  useful  for  thin  fabric  coverings. 
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Test  6:  Three-dimensional  transient  conduction,  constant  heat  flux  impinging  on  a circular  area  (disk) 
on  a semi-infinite  region. 

The  impinging  flux  is  q,  there  is  convective  cooling  from  the  surface,  h(T5  — TJ,  and  the  disk  radius 
is  R.  The  temperature  at  the  center  of  the  heated  area  is  given  by  (Thomas  ref[24],  Eq  5) 


m = ^ 


4at 


K \\  uB} 


1 -exp 


^ f 
h J 


K 


<i)=0 


1 - exp 


-h^R 


4at 

2d2\ 


Ji 


+ erfc 


( R 


4k}uP' 


n 


2 

<0  e " erfc((i)) 


(A12) 


The  first  part  of  the  above  expression  is  the  center  point  temperature  if  there  is  no  convective  cooling. 
Tests  indicate  that  the  accuracy  of  the  FDA  for  this  test  is  primarily  dependent  on  how  accurately  the 
circular  flux  pattern  is  represented  on  the  rectangular  surface  grid.  A small  grid  and  assigning  cell  heat 
gain  according  to  the  portion  of  the  cell  that  is  within  the  circle  improve  accuracy. 

Note  that  tests  1 through  6 involve  a step  change,  which  should  be  the  worst  condition  to  simulate  with 
the  FDA.  In  all  cases  the  maximum  errors  occurred  at  the  first  time  step,  and  the  error  declined  as  time 
increased. 


Test  7:  Three-dimensional  transient  conduction,  uniformly  moving  point-source  heat  flux. 

This  test  consists  of  a point  source  of  power  Q moving  in  an  infinite  body  at  a constant  velocity  v in  the 
X direction  (Schneider,  refI25],  pp.3-86,  Eq  78) 


Kr  1 

— (T-rj  = — exp 
Q " 471  ^ 


2a 


(A13) 


where  | = x - vt  and  t = Adjusting  for  a semi-infinite  body  with  the  point  source 

moving  along  an  adiabatic  surface  is  done  by  replacing  47r  by  27r  in  (A  13).  Generally  good  agreement 
was  achieved  for  this  test.  Accuracy  was  limited  by  grid  size  near  the  point  source  and  the  fact  the  this 
is  a quasi-steady  case  in  that  movement  of  the  point  source  does  not  have  a beginning  point. 


Larkin’s  method 

The  FDA  algorithm  in  CTEST3  was  transferred  directly  into  the  new  substrate  model  TMPSUB2.  The 
addition  of  pyrolysis  forced  the  use  of  a very  small  grid  in  the  region  of  peak  temperature  for  a 
satisfactory  solution.  This  combined  with  the  stability  requirement  of  the  explicit  Euler  method  forced 
a very  small  time  step  and  therefore  a very  long  execution  time.  A different  FDA  algorithm  had  to  be 
found  to  achieve  a program  fast  enough  to  be  useful.  Larkin’s  method  was  chosen  because  of  its 
simplicity  in  that  it  uses  the  same  spatial  FDA  as  the  original  method  while  the  new  temporal  FDA  does 
not  require  the  solution  of  simultaneous  equations. 
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CTEST3  was  not  rewritten  to  run  all  of  the  test  cases,  but  several  tests  were  made  with  TMPSUB2  (some 
using  a modified  surface  boundary  condition)  which  indicate  the  accuracy  of  the  method  for  different  grid 
size  and  time  step  options.  The  results  of  these  tests  are  shown  in  tables  Al,  A2,  and  A3. 

Table  Al  gives  the  results  of  several  tests  which  can  be  compared  to  Eq.(A5)  to  determine  the  effect  of 
grid  spacing.  These  tests  use  a region  40  mm  thick  to  simulate  a semi-infinite  body  which  is  shown  to 
be  appropriate  by  having  negligible  heat  flux  at  the  constant  temperature  surface  at  Z=40  mm. 
Comparisons  are  made  based  on  the  temperature  of  the  surface  (Z=0  mm).  Test  2a:  using  a constant 
1 mm  grid  spacing  the  temperature  after  100  seconds  (502.543)  is  0.33%  less  than  the  exact  value. 
Test  2b:  using  a constant  0.5  mm  spacing  gives  a surface  temperature  0.18%  below  the  theoretical  value. 
Tests  2c  through  2f  use  variable  grid  spacing  to  reduce  the  number  of  cells  and  execution  time  at  the  cost 
of  some  loss  of  accuracy. 

Table  A2  gives  the  results  of  several  tests  where  the  parameters  that  control  the  time  step  while 
maintaining  a constant  grid  spacing,  are  varied.  These  parameters  are  the  maximum  time  step,  dtmax, 
and  the  maximum  temperature  change,  dTmax.  (Whenever  T^+j  - T„  exceeds  dTmax,  the  time  step  is 
halved.)  Obviously  the  greatest  accuracy  should  be  achieved  with  small  values  for  these  two  parameters, 
but  execution  time  is  reduced  by  using  large  values.  There  is  no  obvious  optimum;  the  user  must  choose 
values  appropriate  for  results  he  wishes  to  achieve. 

Table  A3  shows  tests  of  different  grid  spacings  for  full  three-dimensional  heat  conduction  from  a 
stationary  spot  heat  flux.  These  tests  use  representative  thermal  properties  for  the  fabric  and  padding. 
Various  combinations  of  cell  spacings  are  used  to  select  the  best  grid. 
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Table  Al;  Grid  Spacing  Tests 


Thermal  diffusivity:  2e-07  m^/ s 
Thermal  conductivity;  0.1  W/mK 
Surface  heat  flux:  le+04  W/m^ 
Distance  from  surface:  0 m 


time 

Texact 

Test2a 

Test2b 

Test2c 

Test2d 

Test2e 

Test2f 

0.000 

0.0000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.125 

17.8412 

5.000 

9.319 

9.319 

9.319 

9.319 

9.319 

0.250 

25.2313 

9.643 

16.889 

16.889 

16.889 

16.889 

16.889 

0.375 

30.9019 

13.970 

23.331 

23.331 

23.331 

23.331 

23.331 

0.500 

35.6825 

18.298 

28.916 

28.916 

28.916 

28.916 

28.916 

0.625 

39.8942 

22.252 

33.840 

33.840 

33.840 

33.840 

33.840 

0.750 

43.7019 

26.206 

38.035 

38.035 

38.035 

38.034 

38.034 

0.875 

47.2035 

29.838 

42.229 

42.229 

42.229 

42.229 

42.229 

1.000 

50.4627 

33.469 

45.757 

45.757 

45.757 

45.757 

45.757 

1.250 

56.4190 

40.170 

52.366 

52.366 

52.366 

52.366 

52.366 

1.500 

61.8039 

46.383 

58.208 

58.208 

58.208 

58.208 

58.208 

1.750 

66.7558 

52.167 

63.496 

63.496 

63.496 

63.496 

63.496 

2.000 

71.3650 

57.577 

68.361 

68.361 

68.361 

68.361 

68.361 

2.250 

75.6940 

62.656 

73.042 

73.042 

73.041 

73.041 

73.041 

2.500 

79.7885 

67.445 

77.173 

77.173 

77.173 

77.173 

77.171 

2.750 

83.6828 

71.859 

81.305 

81.305 

81.304 

81.304 

81.302 

3.000 

87.4039 

76.273 

85.049 

85.049 

85.048 

85.047 

85.044 

3.500 

94.4070 

84.279 

92.243 

92.243 

92.242 

92.240 

92.233 

4.000 

100.9253 

91.625 

98.911 

98.909 

98.908 

98.904 

98.893 

4.500 

107.0474 

98.433 

105.154 

105.152 

105.149 

105.143 

105.125 

5.000 

112.8379 

104.796 

111.046 

111.042 

111.038 

111.030 

111.002 

6.000 

123.6077 

116.456 

121.976 

121.970 

121.962 

121.946 

121.894 

7.000 

133.5116 

127.018 

132.004 

131.993 

131.981 

132.050 

131.966 

8.000 

142.7299 

136.650 

141.315 

141.300 

141.283 

141.309 

141.188 

9.000 

151.3880 

145.711 

150.034 

150.014 

149.991 

149.993 

149.830 

10.000 

159.5769 

154.229 

158.272 

158.248 

158.219 

158.201 

157.994 

15.000 

195.4410 

191.150 

194.315 

194.266 

194.208 

194.109 

193.654 

20.000 

225.6758 

221.986 

224.670 

224.600 

224.514 

224.343 

223.631 

25.000 

252.3132 

249.025 

251.396 

251.307 

251.195 

250.960 

249.998 

30.000 

276.3953 

273.438 

275.472 

275.367 

275.235 

274.852 

273.660 

35.000 

298.5410 

295.758 

297.599 

297.487 

297.339 

296.828 

295.425 

40.000 

319.1538 

316.566 

318.155 

318.040 

317.879 

317.378 

315.771 

45.000 

338.5137 

336.042 

337.525 

337.408 

337.234 

336.650 

334.844 

50.000 

356.8248 

354.491 

355.815 

355.694 

355.509 

354.928 

352.929 

60.000 

390.8820 

388.738 

389.887 

389.758 

389.550 

388.899 

386.532 

70.000 

422.2008 

420.206 

421.230 

421.093 

420.862 

420.147 

417.431 

80.000 

451.3517 

449.479 

450.409 

450.262 

450.010 

449.235 

446.188 

90.000 

478.7307 

476.960 

477.816 

477.660 

477.388 

476.556 

473.195 

100.000 

504.6265 

502.943 

503.740 

503.574 

503.282 

502.397 

498.738 

cells : 

41 

81 

35 

25 

18 

12 

steps; 

153 

154 

154 

154 

153 

153 

TEST2A:  dz= 

1.0mm,  nz= 

41,  nc=41. 

r=1.000. 

dtmax=l 

. 0 , dTmax 

= 5.0 

TEST2B;  dz= 

0.5mm,  nz= 

81,  nc=81. 

r=1.000. 

dtmax=l 

. 0 , dTmax 

= 5.0 

TEST2C;  dz= 

0.5mm,  nz= 

35,  nc=4. 

r=1.055. 

dtmax=l 

. 0 , dTmax 

= 5.0 

TEST2D;  dz= 

0.5mm,  nz= 

25,  nc=4. 

r=1.115. 

dtmax=l 

. 0 , dTmax 

= 5.0 

TEST2E:  dz= 

0.5mm,  nz= 

18,  nc=4. 

r=1.234. 

dtmax=l 

. 0 , dTmax 

= 5.0 

TEST2F:  dz= 

0.5mm,  nz= 

12,  nc=4. 

r=1.628. 

dtmax=l 

. 0 , dTmax 

= 5.0 
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Table  A2:  Time  Step  Control  Tests 

Thermal  dif fusivity : 2e-07  m^/s 
Thermal  conductivity:  0.1  W/mK 
Surface  heat  flux:  le+04  W/m" 
Distance  from  surface:  0 m 


time 

Texact 

Test2b 

Test2g 

Test2h 

Test2J 

0.000 

0.0000 

0.000 

0.000 

0.000 

0.000 

0.125 

17.8412 

9.319 

9.093 

0.250 

25.2313 

16.889 

16.685 

20.000 

20.000 

0.375 

30.9019 

23.331 

23.148 

0.500 

35.6825 

28.916 

28.750 

31.538 

31.538 

0.625 

39.8942 

33.840 

33.688 

0.750 

43.7019 

38.035 

38.107 

39.690 

39.690 

0.875 

47.2035 

42.229 

42.125 

1.000 

50.4627 

45.757 

45.801 

47.841 

47.841 

1.250 

56.4190 

52.366 

52.390 

1.500 

61.8039 

58.208 

58.223 

58.553 

58.553 

1.750 

66.7558 

63.496 

63.506 

2.000 

71.3650 

68.361 

68.368 

69.264 

69.264 

2.250 

75.6940 

73.042 

72.899 

2.500 

79.7885 

77.173 

77.158 

77.403 

77.403 

2.750 

83.6828 

81.305 

81.191 

3.000 

87.4039 

85.049 

85.028 

85.542 

85.542 

3.500 

94.4070 

92.243 

92.222 

4.000 

100.9253 

98.911 

98.891 

97.831 

97.831 

4.500 

107.0474 

105.154 

105.136 

5.000 

112.8379 

111.046 

111.029 

110.120 

110.120 

6.000 

123.6077 

121.976 

121.963 

120.390 

120.390 

7.000 

133.5116 

132.004 

131.993 

130.660 

130.660 

8.000 

142.7299 

141.315 

141.312 

139.622 

139.622 

9.000 

151.3880 

150.034 

150.053 

148.583 

10.000 

159.5769 

158.272 

158.312 

156.621 

157.545 

15.000 

195.4410 

194.315 

194.402 

192.945 

20.000 

225.6758 

224.670 

224.773 

223.337 

221.615 

25.000 

252.3132 

251.396 

251.504 

250.229 

30.000 

276.3953 

275.472 

275.656 

274.420 

35.000 

298.5410 

297.599 

297.855 

296.726 

40.000 

319.1538 

318.155 

318.512 

317.415 

313.420 

45.000 

338.5137 

337.525 

337.900 

336.887 

50.000 

356.8248 

355.815 

356.236 

355.255 

60.000 

390.8820 

389.887 

390.336 

389.440 

382.464 

70.000 

422.2008 

421.230 

421.690 

420.860 

80.000 

451.3517 

450.409 

450.870 

450.093 

442.530 

90.000 

478.7307 

477.816 

478.274 

477.541 

100.000 

504.6265 

503.740 

504.191 

503.495 

495.723 

cells : 

81 

81 

81 

81 

steps : 

154 

724 

106 

42 

TEST2B: 

dz=0. 5mm, 

nz=81. 

nc=81. 

r=1.000. 

dtmax=1.0. 

dTmax 

= 5.0 

TEST2G: 

dz=0 . 5mm, 

nz=81. 

nc=81. 

r=1.000. 

dtmax=1.0. 

dTmax 

= 1.0 

TEST2H: 

dz=0 . 5mm, 

nz=81 , 

nc=81. 

r=1.000. 

dtmax=1.0. 

dTmax 

= 20.0 

TEST2I: 

dz=0 . 5mm, 

nz=81. 

nc=81 , 

r=1.000. 

dtmax=4. 0, 

dTmax 

= 5.0 

( Same  results  as  test2b  because 

dT  > 2.5  at  dt  = 

1.0  ) 

TEST2J: 

dz=0 . 5mm, 

nz=81. 

nc=81. 

r=1.000. 

dtmax=4 . 0 , 

dTmax 

= 20.0 
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Table  A3 : 

3-D  Transient  Conduction  Tests 

time 

test3a 

test3b 

test3c 

test3d 

test3e 

0.00 

0.000 

0.000 

0.000 

0.000 

0.000 

1.00 

65.749 

70.904 

71.046 

71.254 

71.301 

2.00 

88.450 

91.190 

91.381 

91.660 

91.723 

3.00 

104.807 

106.828 

107.053 

107.382 

107.455 

4.00 

118.536 

120.437 

120.692 

121.060 

121.140 

5.00 

130.745 

132.693 

132.974 

133.375 

133.460 

10.00 

178.269 

180.414 

180.785 

181.287 

181.377 

15.00 

211.444 

213.792 

214.214 

214.765 

214.850 

20.00 

236.140 

238.668 

239.121 

239.692 

239.766 

30.00 

269.762 

272.068 

272.539 

273.119 

273.173 

40.00 

291.304 

293.512 

293.988 

294.550 

294.585 

50.00 

305.937 

308.182 

308.659 

309.201 

309.219 

60.00 

316.279 

318.588 

319.064 

319.589 

319.594 

70.00 

323.830 

326.195 

326.670 

327.180 

327.176 

80.00 

329.498 

331.903 

332.376 

332.876 

332.865 

90.00 

333.853 

336.285 

336.757 

337.248 

337.233 

100.00 

337.271 

339.718 

340.188 

340.673 

340.654 

cells : 

18081 

18081 

18081 

19176 

18375 

steps : 

152 

154 

154 

154 

154 

time: 

303.08 

316.37 

316.10 

338.90 

346.32 

TEST3A: 

dx=0.5mm,  nx=41,  ny=21,  nz  = 21, 
rx=ry=l.i67,  rz=1.187,  dtmax=1.0. 

nc  = 4 , 
dTmax  = 

yw=zw  = 40, 
5.0 

TEST3B; 

dx=0.25mm,  nx=41,  ny=21,  nz  = 21, 
rx=ry=1.240,  rz=1.240,  dtmax=1.0. 

nc  = 4, 
dTmax  = 

yw=zw  = 40 
5.0 

TEST3C: 

dx=0.25mm,  nx=41,  ny=21,  nz  = 21, 
rx=ry=1.210,  rz=1.210,  dtmax=1.0. 

nc  = 4, 
dTmax  = 

yw=zw  = 30 
5.0 

TEST3D: 

dx=0.25mm,  nx=47,  ny=24,  nz  = 17, 
rx=ry=1.161,  rz=1.326,  dtmax=1.0. 

nc  = 4, 
dTmax  = 

yw=zw  = 30 
5.0 

TEST3E: 

dx=0.25mm,  nx=49,  ny=25,  nz  = 15, 
rx=ry=1.149,  rz=1.431,  dtmax=1.0. 

nc  = 4, 
dTmax  = 

yw=zw  = 30 
5.0 
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APPENDIX  B 


TMPSIJB2  User’s  Guide 


General  Information 

TMPSUB2  and  two  helper  programs  are 
available  on  a diskette  for  IBM  PC 
compatible  computers.  The  distribution 
diskette  includes  both  executable  and 
source  code  for  TMPSUB2,  executable 
code  for  the  other  two  programs,  and 
sample  input  files.  The  general 
relationship  of  programs  and  files  is 
illustrated  in  figure  Bl.  The  TSDATA 
program  is  used  to  prepare  data  files  for 
TMPSUB2  which  in  turn  creates  two 
types  of  output  files.  The  plot  file  is  used 
by  the  TSPLOT  program  to  display 
contour  plots  of  substrate  temperatures. 

The  list  file  includes  a step-by-step  record 
of  the  highest  temperature  in  the 
substrate. 

TSDATA  and  TSPLOT  must  be  run  using  MS-DOS  on  IBM  PC  compatible  computers  with  VGA 
graphics.  TMPSUB2  requires  a 386  class  PC  with  math  coprocessor  or  486  class  PC  (no  graphics 
needed)  in  order  to  achieve  satisfactory  performance.  Typical  execution  times  are  20  to  30  minutes  on 
a 33Mhz  486  computer. 

The  TMPSUB2  source  code  (file  TMPSUB2.CCC)  can  be  compiled  using  any  ANSI  C compiler.  All 
TMPSUB2  input  and  output  files  are  ASCII  files.  Therefore,  TMPSUB2  can  be  recompiled  and  run  on 
a different  computer,  while  still  using  TSDATA  and  TSPLOT  on  a PC  and  transferring  files  between  the 
computers.  A different  computer  may  allow  TMPSUB2  to  execute  faster  and/or  handle  more  cells  for 
more  accurate  simulation. 

Be  sure  to  inspect  the  README  file  on  the  distribution  diskette.  One  way  to  read  this  file  is  to  place 
the  diskette  in  drive  A:  (or  drive  B:)  and  type  MORE  < AiREADME  (or  MORE  <B:README).  A 
permanent  copy  may  be  made  with  PRINT  < AiREADME.  The  README  file  contains  a list  of  all  files 
on  the  diskette,  instructions  for  installing  the  necessary  files  on  your  hard  disk,  and  information  on  any 
changes  or  additions  to  the  program. 

There  should  be  at  least  1,000,000  bytes  available  on  your  hard  disk.  It  is  best  to  create  a single 
subdirectory  for  the  executable  programs  and  related  data  files.  This  will  allow  you  to  easily  delete  all 
the  files  related  to  this  program  when  you  are  finished  with  it.  In  general,  when  running  on  a PC,  keep 
all  files  in  the  current  working  direaory. 

The  following  sections  give  details  of  the  operation  of  TMPSUB2,  TSDATA,  and  TSPLOT. 


Figure  Bl.  TMPSUB2  Programs  and  Files 
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TMPSUB2 


All  input  for  TMPSUB2  must  be  placed  in  a data  file.  This  file  is  read  as  the  "standard  input  stream", 
so  TMPSUB2  executes  on  MS-DOS  and  UNIX  computers  by  redirecting  the  input  file.  For  example, 
if  there  already  exists  an  input  data  file  — say,  TEST1.DAT  — one  can  begin  a run  by  typing 
TMPSUB2  < TEST1.DAT  Note:  TMPSUB2  can  be  aborted  at  any  time  by  simultaneously  pressing 
CTRL  and  C. 

Sample  run 

If  you  have  not  already  done  so,  install  the  TMPSUB2  program  following  the  instructions  in  the 
README  file.  For  example,  from  the  directory  on  the  hard  disk  where  you  want  TMPSUB2  installed 
and  with  the  diskette  in  drive  A:,  type  AiINSTALL  A:  . At  the  present  time  you  only  need  to  install  in 
response  to  the  first  question  you  will  be  asked. 

During  installation  a program  will  determine  the  "initial  unallocated  memory".  If  the  required  memory 
is  greater  than  this  amount,  follow  the  instructions  under  memory  requirements  below  before  proceeding. 

Should  you  wish  to  examine  the  directory  at  this  point,  a DIR  command  should  indicate  the  presence  of 
at  least  the  following  files:  TMPSUB2.EXE,  TEST1.DAT,  TEST2.DAT,  and  MEMREM.EXE.  As 
mentioned  above,  begin  the  sample  run  by  typing 
TMPSUB2  < TEST1.DAT 

TMPSUB2  will  then  (1)  indicate  that  it  is  reading  the  data  file,  (2)  echo  the  simulation  title  to  the  screen, 
(3)  initialize  the  data  in  all  cells,  and  (4)  report  the  amount  of  unallocated  memory.  If  there  is 
"insufficient  memory",  an  error  message  will  be  displayed  and  the  run  aborted.  See  the  memory 
requirements  section  below  for  corrective  action. 

Every  time  TMPSUB2  completes  a time  step  during  the  simulation,  it  displays  the  time  (in  seconds)  and 
peak  temperature  (degrees  C)  on  the  screen.  This  allows  you  to  monitor  the  progress  of  the  simulation. 

This  same  information  is  included  on  the  list  file,  TESTl  .LST,  thus  providing  a permanent  record  of  the 
primary  value  computed  during  the  simulation.  The  list  file  also  includes  an  echo  of  the  input  file  which 
helps  to  identify  the  simulation  and  is  especially  useful  in  identifying  any  errors  in  the  input  which  are 
also  written  to  the  list  file. 

The  plot  file,  TESTl.PLT,  contains  all  temperatures  in  the  Y=0  and  Z=0  planes  at  the  times  specified 
plus  data  to  identify  the  simulation  and  the  coordinate  values. 


Memory  requirements 

TMPSUB2  has  been  written  to  handle  an  arbitrary  number  of  cells  up  to  some  limit  imposed  by  available 
memory  or  by  the  operating  system.  MS-DOS  limits  available  random  access  memory  (RAM,  not  disk 
memory)  to  640,000  bytes  which  is  shared  by  the  program,  the  data  in  the  program,  a portion  of  the 
operating  system,  and  perhaps  various  TSR  (Terminate  & Stay  Resident)  programs.  TMPSUB  allocates 
memory  for  its  data  arrays  from  available  RAM.  When  there  is  insufficient  memory  to  run  a particular 
simulation,  you  must  either  make  more  memory  available  or  reduce  the  memory  required  by  the 
simulation.  One  way  to  make  more  memory  available  is  to  remove  TSR  programs,  such  as  network 
connections.  This  usually  requires  rebooting  the  computer  in  such  a manner  that  these  programs  are  not 
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automatically  accessed.  MS-DOS  5.0  uses  somewhat  less  memory  than  previous  versions.  See  you  PC 
consultant  for  assistance.  Run  MEMREM  to  determine  approximately  the  memory  available  for  data. 

The  memory  required  by  the  simulation  is  determined  primarily  by  the  total  number  of  cells  which  in  turn 
is  determined  by  the  values  on  line  4 of  the  data  file  (see  below).  The  total  number  of  cells  can  be 
reduced  to  an  arbitrarily  small  value,  but  accuracy  will  suffer.  If  you  must  search  for  the  number  of  cells 
than  can  be  run  on  your  machine,  begin  with  a relatively  small  number  and  increase  it  until  the 
unallocated  memory  reported  by  TMPSUB2  is  small. 

With  minimal  other  uses  of  memory,  TMPSUB2  is  limited  to  less  than  20,0(X)  cells  under  MS-DOS.  The 
only  way  to  simulate  more  cells  is  to  recompile  TMPSUB2  for  a different  operating  system  and/or  a 
different  computer. 


Contents  of  TMPSUB2  Data  File 


line 

1 

2 

3 

4 

5 


6 


7 


8 

9 


10 


variables  brief  description 

LIST  Name  of  the  output  file 


PLOT 


Name  of  the  plot  file 


TITLE  Project  title;  echoed  to  output 

NX  NY  NZ  number  of  cells  in  the  X,  Y and  Z directions 
( total  number  of  cells  = NX  x NY  x NZ  ) 


NCX  DXO  XW  XO  X-coordinate  data*: 

NCX  number  of  constant  length  nodes 
in  both  directions  from  XO 
DXO  length  [mm]  of  constant  length  nodes 
XW  total  length  [mm]  of  X axis 

XO  location  of  center  of  region  of  constant  length  nodes 
= initial  position  of  peak  flux  from  cigarette 

NCY  DYO  YW  Y-coordinate  data*; 

NCY  number  of  constant  width  nodes  from  Y = 0 
DYO  width  [mm]  of  constant  width  nodes 
YW  total  width  [mm]  of  Y axis 

NCZ  DZO  ZW  FT  Z-coordinate  data*: 

NCZ  number  of  constant  depth  nodes  from  Z = 0 
DZO  width  [mm]  of  constant  depth  nodes 
ZW  total  depth  [mm]  of  Z axis 

FT  fabric  thickness  [mm]^ 

Width  of  air  gap  between  fabric  and  padding  [mm] 

Fabric  data: 

E emittance 

TM  maximum  temperature^  [“C] 

NCF  number  of  data  points  for  T,  K,  C (NCF  <=  10) 

Fabric  thermal  properties'*: 

T temperature  [°C] 

K conductivity  [W/mK] 

C specific  heat  [kJ/kgK] 

Repeat  line  10  up  to  NCF  times. 


SEP 

E TM  NCF 

T K C 
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11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 


DV  DC  DA 

A N1  N2  T H 


A N1  N2  T H 


A N1  N2  T H 
E TM  NCP 

T K C 

DV  DC  DA 

A N1  N2  T H 
A N1  N2  T H 
A N1  N2  T H 
BC  TA  TO  OX 


QO  V Y X+  X- 


HC  HG 


Densities: 

DV  virgin  material  [kg/m^] 

DC  char  [kg/m^] 

DA  ash  [kg/m^] 

Fabric  non-oxidative  pyrolysis  data: 

A reaction  rate  coefficient  [1/min] 

N1  fabric  mass  exponent 

N2  oxygen  concentration  exponent  (must  be  zero) 

T activation  temperature  [K] 

H heat  of  pyrolysis  [kJ/kg] 

Fabric  oxidative  pyrolysis  data: 

A reaction  rate  coefficient  [1/min] 

N1  fabric  mass  exponent 

N2  oxygen  concentration  exponent 

T activation  temperature  [K] 

H heat  of  pyrolysis  [kJ/kg] 

Fabric  char  pyrolysis  data  (similar  to  13) 

Padding  data: 

E emittance 

TM  maximum  temperature^  ["C] 

NCP  number  of  data  points  for  T,  K,  C (NCP  <=  10) 

Padding  thermal  properties'*: 

T temperature  [°C] 

K conductivity  [W/mK] 

C specific  heat  [kJ/kgK] 

Repeat  line  16  up  to  NCP  times. 

Densities : 

DV  virgin  material  [kg/m^] 

DC  char  [kg/m^] 

DA  ash  [kg/m^] 

Padding  non-oxidative  pyrolysis  data  (similar  to  12) 

Padding  oxidative  pyrolysis  data  (similar  to  13) 

Padding  char  pyrolysis  data  (similar  to  13) 

Boundary  & initial  conditions: 

BC  0 = adiabatic  outer  boundaries^, 

1 = constant  temperature  outer  boundaries 
TA  ambient  temperature  ["C] 

TO  initial  substrate  temperature  [“C] 

OX  oxygen  concentration  [fraction] 

Moving  radiant  flux  distribution  on  surface; 

QO  peak  flux  [kW/m^] 

V velocity  [mm/min] 

y standard  deviation  in  Y direction  [mm] 

X+  std  dev  in  positive  X direction  [mm] 

X-  std  dev  in  negative  X direction  [mm] 

Heat  transfer  coefficients^ 

HC  for  quiescent  air  [W/m^K] 

HG  for  impinging  air  [W/m^K] 
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I 


24  dt  DT  TT 


25  LIST  I J K 


26  PLOT 

Notes: 

* The  relationship  between  the  number  of  nodes,  number  of  constant  size  nodes,  size  of  constant  size 
nodes  and  total  length  of  an  axis  must  be  such  that  the  variable  width  nodes  are  increasing  in  size. 

^ The  fabric  thickness  and  the  depth  of  the  constant  depth  cells  are  related:  the  fabric/padding  boundary 
must  be  halfway  between  two  cells.  Therefore,  the  fabric  thickness  must  be  0.5,  or  1.5,  or  2.5,  etc., 
times  the  depth  of  the  constant  depth  cells,  or  the  constant  depth  cells  must  be  2,  or  2/3,  or  2/5,  etc. 
times  the  thickness  of  the  fabric.  The  boundary  between  the  fabric  and  padding  must  be  within  the  region 
of  constant  depth  cells. 

^ When  the  temperature  of  any  fabric  (or  padding)  cell  reaches  the  "maximum"  fabric  (or  padding) 
temperature,  the  simulation  is  terminated. 

^ The  way  that  k(T)  and  c(T)  are  input  is  by  entering  the  values  for  each  at  a number  of  temperatures; 
the  program  then  carries  out  a cubic  spline  fit  to  those  points,  to  obtain  the  values  at  any  other 
temperature.  Even  if  k(T)  and  c(T)  are  given  by  explicit  equations,  this  is  still  the  way  that  the  program 
"knows"  the  values. 

^ The  "outer  boundaries"  of  the  substrate  consist  of  the  X=0,  X=XW,  Y = YW,  and  Z=ZW  planes  (see 
items  5,  6,  & 7). 

^ These  coefficients  refer  to  cases  A and  B in  section  2f.  When  using  HC  set  HG  to  zero  and  vice  versa. 

The  contents  of  the  data  file  are  described  further  in  the  description  of  the  TSDATA  program.  It  will 
generally  be  easier  to  check  a data  file  by  processing  it  with  TSDATA,  than  to  compare  it  line  by  line 
with  the  description  given  above. 


simulation  control: 
dt  maximum  time  step  [s] 

DT  maximum  temperature  change  [°C] 

TT  total  simulation  time  [s] 

Debug  reports  (TMPSUB2  compiled  with  DEBUG  defined) : 
LIST  1 = activate  data  dumps  (default  0) 

I J K are  the  indices  of  the  cell  to  be  studied 

Time  at  which  to  write  temperatures  to  the  plot  file; 
repeat  line  26  up  to  10  times. 
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TSDATA 


The  TMPSUB2  data  file  can  be  created  with  any  ASCII  line  editor.  However,  the  contents  of  the  file 
are  quite  cryptic  and  thus  prone  to  error;  far  better  is  to  use  TSDATA.  TSDATA  is  an  interactive 
program  for  creating  data  files.  Because  it  is  interactive,  it  uses  certain  commands  which  restrict  its 
operation  to  IBM  PC  compatible  computers.  It  includes  extensive  checking  of  the  input  data.  TSDATA 
is  especially  useful  for  creating  a data  file  which  is  only  slightly  different  from  another  data  file.  This 
is  useful  in  performing  the  parametric  studies  for  which  TMPSUB2  was  designed. 

Two  special  files  are  used  by  TSDATA.  The  help  file,  TSDATA. HLP,  contains  the  text  of  the 
interactive  help  messages.  Help  is  activated  by  pressing  the  FI  function  key.  If  the  help  file  is  not 
available  in  the  current  working  directory,  no  interactive  help  will  be  available.  The  configuration  file, 
TSDATA. CFG,  sets  the  colors  of  the  display.  The  file  included  on  the  distribution  diskette  assumes  that 
a standard  VGA  monitor  is  being  used.  If  the  configuration  file  is  not  in  the  current  working  directory, 
a set  of  default  colors  will  be  used.  A new  configuration  file  can  be  made  by  using  the  MAKECFGT 
program.  See  the  README  file  for  instructions. 

The  operation  of  TSDATA  is  explained  on  the  following  pages  which  show  the  messages  and  input 
screens  which  will  appear  as  the  program  is  run.  After  reading  through  these  pages,  try  using  TSDATA 
with  one  of  the  sample  data  files.  Begin  the  program  by  typing  TSDATA.  Abort  the  program  by 
pressing  CTRL  and  C. 

Sample  Runs 

At  this  point  the  following  files  should  be  available  in  the  current  working  directory:  TSDATA.EXE, 
TSDATA. HLP,  TSDATA. CFG,  and  TEST1.DAT.  If  they  are  not,  use  the  INSTALL  procedure. 

First  use  TSDATA  to  view  the  contents  of  the  TEST1.DAT  file.  Type  TSDATA.  Press  ENTER  until 
the  main  menu  appears;  press  ENTER  again  to  set  up  the  file  information.  Enter  TEST1.DAT  as  the 
name  of  the  previous  data  file.  Press  ENTER  to  cycle  through  file  names  (respond  Y to  the  warning 
message  about  a duplicate  file).  Press  ESC  to  return  to  the  main  menu.  Press  ENTER  to  view  each 
section  of  data  in  turn;  press  ESC  to  return.  When  you  reach  the  Save  option,  press  CTRL  and  C to 
abort  the  program.  This  process  will  not  have  changed  TEST1.DAT. 

Now  create  a new  data  file,  TEST3.DAT,  similar  to  TEST1.DAT  but  with  a peak  flux  of  20kW/m^. 
Proceed  as  above  until  you  reach  "Name  of  new  data  file:".  At  this  point  move  the  cursor  to  the  "1"  in 
the  file  name  and  press  3 and  then  press  ENTER.  Revise  the  other  two  file  names  and  the  title  accor- 
dingly. Press  ESC  to  return  to  the  main  menu,  and  move  to  the  boundary  conditions  option.  Press 
ENTER  and  move  (using  cursor  keys)  to  the  peak  heat  flux.  Change  the  value  to  20,  press  ENTER, 
press  ESC.  Now  enter  the  save  option,  and  respond  Y to  save  the  file  TEST3.DAT.  Then  enter  the  exit 
option,  respond  Y to  exit,  and  respond  Y to  create  a batch  file.  After  this  you  should  have  exited 
TSDATA.  Now  type  RUN  to  start  TMPSUB2  using  TEST3.DAT  as  input  data. 
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TSDATA  Input  Screens  and  Help  Messages 


I 


Standard  title  page  with  disclaimer  at  start  of  TSDATA  program. 


TSDATA  — interactive  program  to  prepare  data  files  for  TMPSUB2 

Version  1.0 


Developed  at  the  National  Institute  of  Standards  and  Technology. 

Program  author;  George  Walton 

This  program  is  furnished  by  the  government  and  is  accepted  by  any 
recipient  with  the  express  understanding  that  the  United  States 
Government  makes  no  warranty,  expressed  or  implied,  concerning  the 
accuracy,  completeness,  reliability,  usability,  or  suitability  for 
any  particular  purpose  of  the  information  and  data  contained  in 
this  program  or  furnished  in  connection  therewith,  and  the  United 
States  shall  be  under  no  liability  whatsoever  to  any  person  by  reason 
of  any  use  made  thereof.  This  program  belongs  to  the  government. 
Therefore,  the  recipient  further  agrees  not  to  assert  any  proprietary 
rights  therein  or  to  represent  this  program  to  anyone  as  other  than 
a government  program. 


General  description  of  input  process: 


This  program  assists  you  in  preparing  input  data  files  for  the  TMPSUB2 
program.  It  operates  best  by  reading  an  existing  TMPSUB2  data  file 
which  is  then  modified  to  create  a new  data  file.  A sample  data  file 
is  distributed  with  the  program.  It  can  also  be  used  to  enter  data 
from  scratch.  Several  data  files  can  be  created  in  one  TSDATA  session. 

Data  are  processed  interactively  through  a system  of  data  entry  menus. 
Keyboard  input  is  required  from  the  user  whenever  the  cursor  is  in 
a data  entry  field.  A data  entry  field  is  designated  by  a special 
color  as  is  shown  in  the  lower  right  corner  of  this  screen.  This  is 
a standard  pause  allowing  the  user  to  read  the  screen.  Pressing  any 
key  in  response  will  allow  the  program  to  continue. 

While  the  cursor  is  in  a data  entry  field,  it  will  often  be  possible 
to  get  help  by  pressing  the  FI  function  key.  Help  is  intended  to 
give  additional  information  about  the  data.  For  help  to  work,  the 
TSDATA. HLP  file  must  be  in  the  same  directory  as  the  TSDATA.EXE  file. 
Program  execution  may  be  terminated  when  the  cursor  is  in  a data 
entry  field  by  pressing  CTRL  and  C simultaneously. 


68 


General  description  (continued): 


There  will  usually  be  several  data  entry  fields  on  a single  screen. 
The  field  may  be  blank  or  it  may  present  a default  response.  The 
contents  of  a data  entry  field  (even  all  blank)  may  be  edited. 

That  is,  characters  may  be  overwritten,  deleted,  or  inserted. 

It  may  or  may  not  be  necessary  to  make  an  entry  depending  on  context. 

These  fields  are  ordered  from  top  to  bottom  and  left  to  right. 
Pressing  the  tab  key  or  the  down-arrow  key  moves  the  cursor  to  the 
next  field  to  the  right/down;  The  shifted-tab  or  up-arrow  moves 
to  the  previous  field  to  left/up.  Control-home  or  page-up  moves  the 
cursor  to  the  first  field  on  the  screen.  Control-end  or  page-down 
moves  to  the  last  field.  Movement  between  fields  can  be  done  only 
if  an  entry  is  not  required  and  no  other  keys  have  been  pressed. 

Data  entry  begins  in  the  exchange  mode,  i.e.  the  value  of  the  key 
pressed  replaces  the  character  at  the  cursor.  Pressing  the  insert 
key  will  switch  to  the  insert  mode  (and  from  insert  to  exchange). 

The  delete  key  will  remove  the  character  at  the  cursor.  Control-x 
will  clear  the  entire  data  entry  field.  Move  the  cursor  left  with 
the  left-arrow,  control-left-arrow,  or  home  keys.  Move  the  cursor 
right  with  the  right-arrow,  control-right-arrow,  or  end  keys.  When 
the  data  is  satisfactory,  press  the  ENTER  key. 


Various  checks  are  usually  made  on  each  data  entry.  These  may  produce 
a warning  or  error  message  at  the  bottom  of  the  screen.  The  first  line 
of  this  message  indicates  the  nature  of  the  problem.  The  second  line 
indicates  the  severity  of  the  problem,  the  file  and  line  in  the  TSDATA 
source  code  where  the  error  message  originated,  and  whether  some  help 
may  be  available  by  pressing  the  FI  key. 

Most  errors  will  return  you  to  the  data  entry  field  to  correct  the  input. 
Such  errors  include  invalid  characters  or  numeric  values  outside  certain 
situation-dependent  limits. 

Some  problems  may  generate  a question; 

Question?  (y/n)  ^ 

The  user  can  press  the  Y and  then  ENTER  keys  to  indicate  a positive  (yes) 
response  or  N and  then  ENTER  to  indicate  a negative  (no)  response. 

Some  errors  are  fatal  causing  the  program  to  terminate. 
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Screen  1:  (primary  menu) 


TMPSUB2  data  preparation: 

^ File  information 
^ Geometry  description 
^ Fabric  data 
^ Padding  data 
^ Boundary  conditions 
^ Simulation  control 
^ Save  this  data  file 

^ Exit  data  preparation 


Use  cursor  keys  to  move  between  menu  selections. 

Press  ENTER  to  activate  the  menu  selection  at  the  X. 

Press  ESC  to  return  from  a selection.  Press  FI  for  help. 


General  help  message: 

This  program  assists  you  in  preparing  input  data  files  for  the  TMPSUB2 
program. 

This  is  the  main  menu.  It  directs  you  to  the  different  data 
preparation  subsections. 

Begin  each  data  file  by  entering  the  "file  information". 

Data  relating  to  the  geometry,  fabric,  padding,  or  boundary  conditions 
may  be  entered  or  changed  in  any  order. 

The  data  file  is  not  created  until  you  "save  this  data  file". 

Multiple  data  files  can  be  created  in  one  interactive  session. 
Terminate  the  session  by  entering  "exit  data  preparation". 

Note: 

An  X appears  in  the  selection  designator: 
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Screen  2: 


File  information: 


press  ESC  when  done;  press  FI  for  help. 


Copy  data  from  previous  data  file  named: 
Name  of  new  data  file: 

Name  of  list  file: 

Name  of  plot  file: 

Title  for  this  simulation: 


Files  previously  defined  in  this  session: 
data  files  list  files  plot  files  titles  (first  40  characters) 


General  help  message: 

The  data  files  you  create  are  stored  in  the  current  directory.  If  you 
want  to  make  a data  file  which  is  similar  to  a previous  data  file,  enter 
the  name  of  the  previous  file.  That  file  must  also  be  in  the  current 
local  directory. 

TMPSUB2  creates  two  files  when  it  is  run: 

(1)  the  list  file  (which  echoes  the  input,  saves  any  error  messages,  and 
notes  the  peak  temperature  of  the  substrate  during  the  simulation) , 

(2)  the  plot  file  (which  becomes  the  input  file  to  the  plotting  program 
used  to  view  the  temperature  distribution  in  the  substrate) . 

Files  should  generally  have  different  names.  For  example,  consecutive 
TMPSUB2  runs  with  the  same  plot  file  names  will  save  only  the  plot  file 
from  the  last  run  (which  will  have  replaced  the  previous  plot  files). 

Use  identical  names  if  you  definitely  want  to  replace  existing  files, 
including  data  files. 

The  title  that  you  enter  here  is  echoed  in  the  list  and  plot  files.  Use 
the  title  as  a reminder  of  the  special  features  of  a particular  run. 

Up  to  10  data  files  created  previously  in  this  session  are  recorded  at 
the  bottom  of  the  screen  for  your  convenience. 
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Help  messages  in  response  to  error  messages: 

(1)  could  not  open  file: 

The  file  which  you  specified  was  not  found.  A typing  error  is  the  most 
likely  problem,  or  the  file  is  not  in  the  current  directory. 

(2)  duplicate  file  names: 

This  message  indicates  that  a file  name  matches  one  in  the  local  directory, 
or  one  previously  set  in  this  session,  or  that  you  have  not  pressed  ENTER 
at  each  of  the  file  names  to  complete  the  check  of  file  names. 

If  you  have  copied  an  existing  data  file,  the  "new"  file  names  are 
the  ones  that  appeared  in  the  existing  data  file. 

Files  should  generally  have  different  names.  For  example,  consecutive 
TMPSUB2  runs  with  the  same  plot  file  names  will  save  only  the  plot  file 
from  the  last  run  (which  will  have  replaced  the  previous  plot  files). 

Use  identical  names  if  you  definitely  want  to  replace  existing  files, 
including  data  files. 

You  may  use  a duplicate  file  name  by  responding  'Y'  to  the  question 
about  writing  over  the  previous  file. 


Notes: 

An  entry  for  the  previous  data  file  causes  its  data  to  be  copied. 

The  new  data  file  will  be  created  when  "save  this  data  file"  is  executed  at  the  main  menu.  This  save  adds 
the  new  data,  list  and  plot  file  names  to  the  list  displayed  at  the  bottom  of  this  screen. 

There  should  be  no  duplicate  file  names.  That  would  cause  files  to  be  overwritten  during  simulation. 
The  program  checks  for  duplicates. 

The  simulation  title  is  echoed  in  the  output  files. 

Note  that  upon  entering  the  "Title  for  this  simulation, " the  cursor  jumps  back  to  the  first  box  and  blanks 
it  out.  This  is  done  only  so  that  if  the  user  has  changed  his  mind  (or  an  error  has  been  made),  new  (or 
corrected)  entries  can  be  made  immediately.  If  there  has  been  no  error,  then  simply  press  ESC. 
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Screen  3: 


Geometric  description; 
Fabric  thickness: 

Width  of  air  gap; 


mm 

mm 


press  ESC  when  done; 


press  FI  for  help. 


X-coordinate  data; 

Number  of  grid  points; 
Total  length  of  substrate: 
Size  of  constant  length  cells: 
Number  of  constant  length  cells: 

y-coordinate  data; 

Number  of  grid  points; 
Half-width  of  substrate: 
Size  of  constant  width  cells: 
Number  of  constant  width  cells; 

Z-coordinate  data; 

Number  of  grid  points; 
Total  depth  of  substrate; 
Size  of  constant  depth  cells; 
Number  of  constant  depth  cells; 


^ Rx  = l.xxx 

mm 
mm 

^ Ry  = 1 . XXX 

mm 
mm 

^ Rz  = l.xxx 

mm 
mm 

cells:  xxxxx 


General  help  message: 

The  substrate  consists  of  a thin  fabric,  padding,  and  possibly  an 
air  gap  between  them.  An  air  gap  width  of  0.0  indicates  no  gap. 

The  X direction  is  along  the  cigarette;  the  Y direction  is  along  the 
fabric;  the  Z direction  is  into  the  padding. 

A variable  grid  is  used.  It  consists  of  several  constant  width  cells 
near  the  point  of  peak  incident  heat  flux  followed  by  increasingly 
larger  cells  out  to  the  boundaries  of  the  substrate.  For  example; 

I I I I I I I I I I I 

I I I I I I I I I I I 

Y=0  Y=width 

Each  variable  width  cell  is  R times  longer  than  the  preceding  cell. 
Simulation  is  most  accurate  when  R is  only  slightly  greater  than  one. 

A value  of  R less  than  about  1.25  should  be  sufficiently  accurate. 
Simulation  will  generally  involve  a trade-off  of  accuracy  and  run  time. 

The  fabric  thickness  and  the  depth  of  the  constant  depth  cells  are 
related:  the  fabric/padding  boundary  must  be  halfway  between  two  cells. 

Therefore,  the  fabric  thickness  must  be  0.5,  or  1.5,  or  2.5,  etc.,  times 
the  depth  of  the  constant  depth  cells. 


Other  help  messages: 

For  the  X-axis  the  following  conditions  should  apply: 

4 * NCX  < NX  and  (NX  - 1)  * DXO  <=  XW 
where  NX  = number  of  grid  points,  XW  = total  length  of  substrate, 

DXO  = size  of  const  length  cells,  & NCX  = number  of  const  length  cells. 
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I 


For  the  Y-axis  the  following  conditions  should  apply: 

2 * NCY  < NY  and  (NY  - 1)  * DYO  <=  YW 

where  NY  = number  of  grid  points,  2YW  = total  width  of  substrate,  | 

DYO  = size  of  const  width  cells,  & NCY  = number  of  const  width  cells.  ' 

Because  of  the  assumed  bilateral  symmetry  of  the  flux,  it  is  only  necesary 
to  make  calculations  between  y = 0 and  y = YW  (half  the  width).  The  number  i 

of  constant-width  cells  is  that  in  the  half-width  section.  i 

For  the  Z-axis  the  following  conditions  should  apply:  I 

2 * NCZ  < NZ  and  (NZ  - 1)  * DZO  <=  ZW  I 

where  NZ  = number  of  grid  points,  ZW  = total  depth  of  substrate,  j 

DZO  = size  of  const  depth  cells,  & NCZ  = number  of  const  depth  cells.  I 

The  fabric  thickness  and  the  depth  of  the  constant  depth  cells  are  j 

related:  the  fabric/padding  boundary  must  be  halfway  between  two  cells.  j 

Therefore,  the  fabric  thickness  must  be  0.5  or  1.5  or  2.5  etc.,  times  i| 

the  depth  of  the  constant  depth  cells,  or  the  constant  depth  cells  must  ji 

be  2 or  2/3  or  2/5  etc.  times  the  thickness  of  the  fabric.  |l 

The  boundary  between  the  fabric  and  padding  must  be  within  the  region  |i 

of  constant  depth  cells.  I 

P 

Notes:  j 

I 

"Rx",  "Ry",  and  "Rz"  are  the  geometric  progression  rates.  The  sample  problems  have  a relatively  high  j 

value  in  the  Z direction.  Tests  have  indicated  this  is  satisfactory  because  the  low  conductivity  of  the  i 

padding  permits  less  heat  transfer  in  this  direction.  ! 

"cells:"  gives  the  total  number  of  cells  used  to  model  the  substrate.  It  equals  NX  x NY  x NZ.  This  | 

value  is  critical  to  the  total  memory  required  for  a simulation. 

The  following  figure  illustrates  the  substrate  coordinate  system  and  the  variable  grid: 
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Densities:  ( kg/m^  ) 
Virgin  material: 


char : 


ash : 


General  help  message; 

Emittance  is  used  in  computing  radiant  heat  loss  from  the  fabric  to 
ambient  and  radiant  heat  transfer  across  the  air  gap. 

The  simulation  stops  when  the  temperature  of  any  fabric  cell  reaches 
the  prescribed  maximum  temperature. 

Thermal  conductivity  (K)  and  heat  capacity  (C)  can  be  specified  at  up 
to  10  temperatures.  TMPSUB2  uses  a cubic  spline  curve  fit  for  K and  C. 
As  fabric  mass  is  lost  by  pyrolysis,  K is  also  reduced  in  proportion 
to  the  current  density.  Values  are  required  at  least  two  temperatures 
even  if  K and  C are  constant. 

A two-stage,  three-reaction  pyrolysis  model  is  used:  the  virgin 
material  is  converted  to  char  in  stage  1 and  then  to  ash.  The  char 
density  is  the  value  for  completely  converting  the  virgin  material 
to  char  with  no  conversion  to  ash.  The  ash  density  is  what's  left 
after  total  pyrolysis. 
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The  pyrolysis  equations  are  of  the  form; 

R = A * Dm^nl  * Cx^n2  * exp(  Ta/T  ) and  Q = R * He 
where 

R = rate  of  pyrolysis 
A = reaction  coefficient  [l/min] 

Dm  = density  of  material  (virgin  or  char  ) 

nl  = related  exponent 

Cx  = oxygen  concentration 

n2  = related  exponent 

Ta  = activation  temperature  (°K) 

T = current  cell  temperature  [°K] 

He  = heat  of  pyrolysis  [kJ/kg] 

The  first  reaction  is  thermal  degradation  (n2  = 0 and  He  < 0)  of  virgin 
material.  The  second  reaction  is  oxidation  of  the  virgin  material. 
These  two  reactions  produce  char.  The  third  reaction  is  oxidation 
of  the  char  to  produce  ash. 


Help  message  in  response  to  error  message: 


Temperatures  must  be  given  in  increasing  order. 

The  lowest  temperature  should  be  several  degrees  less  than  any 
possible  substrate  temperature.  The  highest  temperature  should 
be  several  degrees  higher  than  the  maximum  temperature. 


Screen  5: 


Padding  Properties; 


press  ESC  when  done;  press  FI  for  help. 


Emittance; 


Maximum  temperature; 


# 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Densities;  ( kg/m^  ) 
Virgin  material; 


char; 


ash; 


Padding  help  messages  are  similar  to  fabric  messages. 
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Screen  6: 


Boundary  Conditions: 

At  X=0,  X=Xmax,  Y=ymax,  Z=Zmax: 


press  ESC  when  done; 

(^)  adiabatic 

(P)  constant  temperature 


Initial  substrate  temperature: 

Ambient  temperature: 
Oxygen  mass  fraction: 


Heat  transfer  coefficients: 

Quiescent  air; 
Impinging  air; 


press  FI  for  help 


Moving  heat  flux  pattern  from  the 
Peak  heat  flux: 
Initial  X position  of  peak; 

+X  velocity: 

±Y  standard  deviation  (A) ; 

+X  standard  deviation  (B); 

-X  standard  deviation  (C) ; 


arette: 


kW/m^ 

mm 

mm/min 

mm 

mm 

mm 


General  help  message: 

Select  the  condition  of  the  outer  boundaries  (X=0,  X=length,  Y=width, 
Z=depth)  by  pressing  ENTER  in  the  appropriate  place. 

The  incident  heat  flux  from  the  cigarette  is  represented  by  a moving 
flux  distribution  on  the  surface  of  the  substrate.  This  distribution 
has  the  shape  of  non-symmetric  Gaussian  curve  which  is  described  by  6 
parameters : 

(1)  the  X coordinate  at  the  peak  of  the  curve  at  the  start  of  the 

simulation  { the  initial  position  of  the  peak  is  (Xo,  0,  0)  }, 

(2)  the  speed,  S,  at  which  the  peak  moves  along  the  X axis  { the 
position  P of  the  peak  at  time  t is  (Xo+S*t,  0,  0)  }, 

(3)  the  maximum  heat  flux  at  the  peak  { at  position  (P,  0,  0)  }, 

(4)  the  width  of  the  curve  in  the  Y direction  { at  position  (P,  ±A,  0), 
the  heat  flux  is  0.37  times  the  flux  at  the  peak  }, 

(5)  the  width  of  the  curve  in  front  of  the  peak  { at  position  (P+B,  0,  0) 
the  heat  flux  is  0.37  times  the  flux  at  the  peak  },  and 

(6)  the  width  of  the  curve  behind  the  peak  { at  position  (P-C,  0,  0), 
the  heat  flux  is  0.37  times  the  flux  at  the  peak  }. 
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Screen  7: 

simulation  Control: 

Maximum  time  step: 

Maximum  temperature  change: 

Total  simulation  time: 

Plot  times  (s); 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


press  ESC  when  done; 


press  FI  for  help. 


General  help  message: 

TMPSUB2  uses  a variable  time  step  which  is  chosen  so  as  not  to  exceed 
the  maximum  time  step  given  and  so  that  the  maximum  temperature  change 
in  any  cell  is  not  greater  than  the  maximum  given.  This  allows  the 
program  to  run  cjuickly  when  the  temperatures  are  not  changing  rapidly. 
Smaller  values  for  these  two  parameters  will  lead  to  more  accurate 
simulations  at  the  cost  of  longer  execution  times. 

The  simulation  will  stop  at  the  total  simulation  time  unless  it  has 
already  stopped  by  exceeding  a maximum  fabric  or  padding  temperature. 

Plot  times  are  entered  in  increasing  order.  When  the  simulation  reaches 
a plot  time,  substrate  temperatures  are  copied  to  the  plot  file.  Leave 
later  positions  blank  if  you  do  not  want  to  use  all  10  times.  A plot  is 
automatically  written  when  the  simulation  stops.  The  last  plot  time 
should  be  greater  than  the  simulation  time  to  satisfy  a requirement  in 
the  TMPSUB2  program. 
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Screen  8: 


Save  this  data  file?  (y/n)  ^ 
display  any  error  inessage(s). 


General  help  message: 

If  you  entered  this  area  accidentally  and  are  not  ready  to  stop 
preparing  the  current  data  file,  respond  'N'  to  this  question. 
This  will  return  you  to  the  main  menu. 


Help  messages  in  response  to  error  message: 

(1)  check  data. 

Data  in  the  indicated  section  may  be  incomplete  or  incorrect. 

Enter  that  section  and  check  the  data. 

(2)  temperature  problem. 

The  lowest  temperature  for  fabric  or  padding  thermal  properties  must 
be  several  degrees  lower  than  either  the  initial  or  ambient  temperature. 
You  need  to  change  (at  least)  one  of  those  values. 

(3)  files  limit. 

No  more  TMPSUB2  data  files  can  be  created  in  this  session. 

You  probably  should  exit  TSDATA  now.  You  may  replace  one  of  the 
data  files  already  created. 
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Screen  9: 


Exit  data  preparation?  (y/n)  ^ (1) 
display  error  messages,  if  applicable  (2) 
prepare  DOS  batch  file  to  run  all  cases  (y/n)?  ^ (3) 


Help  messages: 

(1)  initial  message. 

If  you  entered  this  area  accidentally  and  are  not  ready  to  stop 
preparing  data  files,  answer  'N'  to  this  question. 

(2)  error  messages. 

This  warning  indicates  that  TSDATA  may  have  data  which  you  have  not  saved. 
If  you  respond  'Y'  to  the  question  about  checking  this  data  file,  you 
are  returned  to  the  main  menu  where  you  can  check  and  then  save  the  file. 
If  you  respond  'N',  you  will  continue  to  exit  the  TSDATA  program. 

(3)  batch  file  message. 

A positive  response  will  create  a file  RUN. BAT  resembling: 

TMPSUB2  <datafile.l 
TMPSUB2  <datafile.2 
TMPSUB2  <datafile.3 
TMPSUB2  <datafile.4 

This  file  can  be  used  to  run  TMPSUB2  on  computers  using  MS-DOS. 


Notes: 

Message  (1)  appears  automatically. 

Error  message  (2)  may  appear. 

Message  (3)  appears  after  a Y to  question  (1). 
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TSPLOT 


The  TSPLOT  program  will  display  contour  plots  of  substrate  temperatures  on  the  screen.  It  is  an 
interactive  program  which  uses  the  plot  file  from  TMPSUB2  as  input.  Because  it  is  interactive  and 
graphic,  it  uses  certain  commands  which  restrict  its  operation  to  IBM  PC  compatible  computers.  A plot 
file  contains  the  substrate  surface  (Z  = 0 plane)  and  center-plane  (Y  = 0 plane,  see  fig  B2)  temperatures 
at  various  times  set  in  the  input  data  file. 

Three  special  files  are  used  by  TSPLOT.  The  CHRSET.VGA  file  contains  the  bit  patterns  for  the 
graphic  display.  This  file  is  required.  The  help  file,  TSPLOT. HLP,  contains  the  text  of  the  interactive 
help  messages.  Help  is  activated  by  pressing  the  FI  function  key.  If  the  help  file  is  not  available  in  the 
current  working  directory,  no  interactive  help  will  be  available.  The  configuration  file,  TSPLOT. CFG, 
sets  the  colors  of  the  display.  The  file  included  on  the  distribution  diskette  assumes  that  a standard  VGA 
monitor  is  being  used.  If  the  configuration  is  not  in  the  current  working  directory,  a set  of  default  colors 
will  be  used.  These  colors  may  not  provide  satisfactory  contour  plots.  A new  configuration  file  can  be 
made  by  using  the  MAKECFGG  program.  See  the  README  file  for  instructions. 

The  operation  of  TSPLOT  is  briefly  explained  on  the  following  page.  Because  the  screen  interface  is 
similar  to  TSDATA  and  relatively  few  options  are  available,  a detailed  description  should  not  be 
necessary.  After  reading  this  description,  try  using  TSPLOT  with  the  sample  plot  files.  Begin  the 
program  by  typing  TSPLOT.  Abort  the  program  by  pressing  CTRL  and  C. 

Sample  Run 

At  this  point  the  following  files  should  be  available  in  the  current  working  directory:  TSPLOT.EXE, 
CHRSET.VGA,  TSPLOT.HLP,  TSPLOT.CFG,  TESTl.PLT,  and  TEST2.PLT.  If  they  are  not,  use  the 
INSTALL  procedure.  If  you  ran  TESTl  or  TEST2  and  aborted  before  completion,  use  INSTALL  to 
replace  the  plot  files  with  the  original  complete  versions. 

Type  TSPLOT.  Press  ENTER  until  the  main  menu  appears;  press  ENTER  again  to  get  the  plot  file. 
Enter  TESTl.PLT  as  the  name  of  the  plot  file.  After  returning  to  the  main  menu,  press  ENTER  to 
display  the  next  plot.  Note  that  this  plot  is  for  t = 1 second,  and  the  pattern  is  rather  small.  After 
viewing  this  plot,  press  any  key  to  return  to  the  main  menu.  Move  (use  cursor  keys)  to  the  set  display 
parameters  option  and  press  ENTER.  Change  the  Xmin  value  from  "0"  to  "10",  press  ENTER,  change 
Xmax  from  "60"  to  "50",  press  ENTER,  and  press  ESC.  Now  press  ENTER  to  display  the  prior  plot. 
The  limits  of  the  display  have  been  changed,  and  the  heated  area  appears  larger.  Press  ENTER  to 
alternate  between  the  main  menu  and  the  next  plot.  In  the  later  plots  note  the  discontinuity  in  the  profiles 
in  the  padding  caused  by  the  air  gap.  After  the  last  plot  you  will  arrive  at  the  exit  option.  Instead  of 
exiting,  move  to  the  get  plot  file  option,  press  ENTER,  and  enter  TEST2.PLT  as  the  name  of  the  plot 
file.  Return  to  the  main  menu  and  display  the  next  plot  which  is  the  first  plot  from  TEST2.PLT.  Note 
that  the  limits  of  the  X-axis  have  not  been  reset  for  this  new  plot.  As  you  continue  by  displaying 
successive  plots  note  how  the  point  of  peak  temperature  is  moving  along  the  X-axis.  When  you  reach 
the  exit  option,  press  ENTER,  and  respond  Y to  exit  the  program. 
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TSPLOT  Main  Menu 


This  is  the  main  menu  screen.  It  directs  you  to  the  different  options. 


TSPLOT 


program: 

Get  plot  file 
Display  next  plot 
Redisplay 

Set  display  parameters 
Exit 


Use  cursor  keys  to  move  between  menu  selections. 

Press  ENTER  to  activate  the  menu  selection  at  the  X. 

Press  ESC  to  return  from  a selection.  Press  FI  for  help. 


The  first  option  is  "Get  plot  file".  You  must  get  a (enter  the  name  of  an  existing)  plot  file  before  any 
plots  can  be  displayed. 

The  plot  file  is  read  sequentially.  You  may  display  the  data  at  the  next  plot  time  (option  two),  or  you 
may  redisplay  the  last  plot  shown  (option  three). 

The  graphic  display  consists  of  three  parts:  the  temperature  scale  (in  °C)  to  the  right,  the  substrate 
surface  temperature  contours  in  the  upper  left  of  the  screen,  and  the  center-plane  temperatures  in  the 
lower  left.  In  other  words,  the  temperature  contours  for  the  Z = 0 and  Y = 0 planes  are  displayed  together 
as  if  they  have  been  folded  along  the  X-axis  so  as  to  both  be  flat  on  the  screen.  The  display  also  shows 
the  time  and  the  X-axis  with  coordinates  at  the  left  and  right  edges  and  tic  marks  every  millimeter. 

The  fourth  option  lets  you  reset  the  display  parameters.  You  may  set  the  coordinates  of  the  left  and  right 
edges  of  the  display.  These  values  change  both  the  position  and  the  scale  of  the  regions  being  displayed 
thus  enlarging  or  shrinking  the  plot.  The  initial  limits  of  the  X-axis  are  for  the  entire  region  simulated. 
You  may  also  change  the  temperature  contours  which  are  initially  set  at  25°C  intervals.  Changes  must 
be  made  so  that  each  temperature  is  always  less  than  the  one  above  and  more  than  the  one  below.  You 
will  then  probably  want  to  redisplay  the  last  plot. 

After  displaying  one  set  of  plots,  you  may  get  a new  plot  file  without  exiting  this  program.  This  is  the 
only  way  to  go  back  in  time  on  the  plot  file:  get  the  plot  file  again  and  start  at  the  beginning. 

When  you  reach  the  end  of  the  plot  file,  TSPLOT  automatically  takes  you  to  the  exit  option. 
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APPENDIX  C 


Analysis  of  numerical  errors  produced  by  having  a runaway  reaction  rate 

We  have  made  an  effort  to  make  a fine  grid  in  the  region  where  the  fluxes  and  flux  gradients  are  high. 
Nevertheless,  when  the  reaction  rates  become  very  high,  the  gradients  will  become  so  steep  that  the 
approximation  of  constant  temperature  and  constant  reaction  rate  within  a cell  becomes  questionable.  In 
this  appendix,  we  examine  the  magnitude  of  the  errors  thus  committed  by  discretization.  This  analysis 
can  also  serve  to  make  appropriate  corrections  in  the  program;  this  has  not  been  done  here,  partly  because 
the  analysis  should  first  be  generalized  to  the  non-symmetric  case  (see  the  assumptions  made  just  below 
Eq.(Cl)). 

Consider  the  heat  diffusion  equation,  Eq.(12).  Suppose  we  have  the  (correct)  temperature  distribution 
T(x,y,z,t),  with  a peak  at  the  reaction  rate  is  given  by  an  expression  such  as  (20)  or  (22).  Let 

us  simplify  this  form  and  assume  that 

Rp  = Kcxp(-TJT(r))  (Cl) 

where  all  the  preexponential  factors  are  lumped  together  as  the  factor  and  r s (x,y,z)  is  the  position 
vector.  Consider  the  terms  in  Eq.(12);  we  simplify  the  analysis  by  assuming  that  the  temperature  peak 
lies  at  the  cell  center  (i.e.,  at  Pq),  and  that  the  distribution  is  symmetric  fore-and-aft.  Then,  since 
calculating  Rp  in  that  cell  means  calculating  it  at  the  center,  and  therefore  the  peak,  that  means  we 
overestimate  the  reaction  rate,  since  the  rate  falls  off  at  the  faces  of  the  cell. 

Symmetry  implies  that 


V-(KVr)  = K 


' d^T  ^ d'^T^ 

dy^  dz^ 


= 3k 


d^T 

dx^ 


A second-order  approximation  to  this  derivative  is 


V(KVr)  = 3k 


(Ax)2 


Because  of  the  assumed  symmetry,  Ti_,  = Tj^.!,  so  that  we  finally  have 


V(KVr)  = 6k 


T -T 

*i-l 

(Ax)^ 


Inserting  this  into  Eq.(12),  that  equation  becomes 

pc 


T*  T 

— = HR(r)  - 2k  - ... 


(C2) 


(C3) 


(C4) 


(C5) 


minus  similar  terms  in  Ay  and  Az 

(as  indicated  by  the  ellipsis).  Suppose  further,  for  the  sake  of  simplicity,  that  the  temperature  profile  is 
Gaussian: 


where 


0(x,y,x,O  = e„exp 


0 ^ T - T, 


(C6) 
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and  = ambient  (reference)  temperature. 

If  we  redefine  the  origin  to  be  at  (Xo,yo,Zo),  then 

e(x,y,x,t)  = e^expC-r^/o^) 

where  r^  = + z^. 

Then  ]f  x^/tr^  < < 1 everywhere  within  the  cell,  we  can  expand  (C7):  at  the  center  of  the  cell, 


(C7) 


01  = e('-o)  = 00^  = 


and 


0.-1  = exp  [-(Ax/ 0)2]  = e„[l-(Ax/o)"]. 


Define  5 = (Ax/o)^  . Then 


so  that  Eq.(C5)  becomes 


7’,  - 7’,.!  = e,  - = 0.^ 


pc—  = hr  - ft 

Af  ^ (Ax)2 

Assuming  three-fold  symmetry,  the  other  two  terms  in  (Cll)  are  the  same,  and  thus 

pc-^  = H^R  -6k  — 


(C8) 


(C9) 


(CIO) 


(Cll) 


(C12) 


Note  that  the  dependence  on  Ax  has  dropped  out.  Indeed,  the  last  term  is  the  exact  diffusive  loss  rate, 
resulting  from  the  distribution  (C7). 

Now  let  us  examine  the  magnitude  of  the  error  made  by  using  a finite  value  of  Ax:  from  Eq.(C7), 

AS  . 6,-6,.,  = e„[l-exp(-?)]  (C13) 

Then  expanding  to  the  next  higher  order  than  was  done  in  (C9), 

Ae=e^ai-5/2)  (C14) 

and  the  second  term  in  (C5)  becomes 

2KAe  2k0„ 


(1-5/2) 


(C15) 


(Ax)2  a2 

Thus  the  numerical  calculation  underestimates  the  loss  rate.  Since  the  correct  heat  loss  rate  is 

= 6K6Ja^  (C16) 

and  we  have 

6k0. 


L = 


(Ax)^ 


[l-exp(0] , 


(C17) 
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we  evidently  must  multiply  L by  the  ratio 


[l-exp(-5)]  1-5/2 

in  order  to  get  the  correct  loss  rate  5£. 


(C18) 


As  the  temperature  "runs  away,"  the  distribution  becomes  increasingly  peaked,  and  o declines.  As  it  does 
so,  ^ increases,  and  so,  therefore,  does  F;  this  helps  to  slow  down,  or  moderate,  the  runaway.  ^ is 
readily  found  to  be 


5 = in(e,/e,.,)  (C21) 

Next,  consider  the  source  term  in  Eqs  (C5)  and  (C12),  again.  Define 

f(.x,y^)  = exp[-r^/r(x,y,z)]  (C22) 

The  correct  power  output  in  the  elemental  volume  AV  is,  from  Eq.  (Cl), 

Ajc/2  Ay/2  A2/2 

RpAV  = j dx  I dy  j f(x,y^)dz  (^21) 

0 0 0 

where  we  have  again  translated  (Xo,yo,Zo)  to  be  at  the  origin,  and  have  taken  advantage  of  the  3-fold 
reflection  symmetry.  Again  assume  Eq.(C7)  to  hold.  If  ^ < < 1,  then  within  the  cell  we  can  Taylor 
expand  and  write 


f(x,yyZ)  = 


dx 


.2fg2f\ 


/O 


/o 


Since  the  peak  is  at  the  origin. 


and  dropping  the  higher-order  terms. 


plus  similar  terms  in  y and  z. 


dx 


df] 


dy 


/o 


. i|i . 


(C22) 


(C23) 


i! 

f 

2 

'ay 

faV] 

2 

iax^J 

p 

0 2 

(C24) 


We  find  that 


Thus 


W)o 


IT,  0 


^exp(-Vrj 


fix,y,z)  = exp(-r^/r^) 


1 

Tja^ 


and  integration  over  the  cell  yields 


(C25) 


(C26) 
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(C27) 


RAV  = (Ax)3/?,/(0) 


1 — ^ c2 

472 


This  is  only  valid  for  Ax/a  < < 1.  In  that  limit,  however,  (C26)  is  well  approximated  by 

Rp  = -R.exp(-r^/r„)exp 


r2 

472 


(C28) 


which  is  a reasonable  first  approximation  to  the  correct  integral,  even  when  Ax/a  is  not  very  small! 
The  numerical  approximation  is 


^p  = ^p(^o)  = ^-exp(-r^/r„) 


(C29) 


Thus  (C28)  shows  that  we  must  multiply  this  approximation  by 
better  approximation  to  the  source  term. 


exp 


^A^jn  ^2 


47 


in  order  to  get  a still 
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APPENDIX  D 


Analysis  of  Ignition  Experiment 

Among  the  things  we  wish  to  extract  from  this  experiment,  is  the  heat  transfer  coefficient  to  the  surface, 
h.  Fig.  16  shows  the  radiative  and  total  heating  fluxes  to  the  gauge,  as  a function  of  the  power  input. 
Thus,  consider  the  18  kW/m^  case.  Without  purge  flow,  0^^*  = 10  kW/m^;  with  the  purge  flow  on, 
= 18  kW/m^,  suggesting  that  the  convective  flux  is  0<,on  ” ^tot  “ ‘^rad  “ 18  - 10  = 8 kW/m-. 
However,  when  the  purge  flow  is  turned  on,  it  cools  off  the  heating  element  somewhat.  How  much  it 
does  so,  however,  is  unknown.  Assume  that  the  purge  flow  reduces  the  hot-surface  temperature  such  that 
(^rad  is  reduced  to  some  fraction  F of  the  original.  Then 

4.^  = lOF  (Dl) 

and 

= 18  - lOF  (D2) 

At  the  point  P,  on  the  heater  axis  but  at  the  surface  being  heated  (5.4  mm  below),  the  heater  subtends 
a solid  angle  such  that  the  view  factor  is  fi.  Hence  the  impinging  radiative  flux  is 

<l>™a  = (1-0)07/  (D3) 

where  is  the  emissivity  of  the  device  surface.  For  the  particular  case  that  was  carried  out,  the  standoff 
distance  is  5.4  mm,  while  the  diameter  of  the  glowing  filament  is  13  mm.  According  to  Siegel  and 
Howell  (ref  [25],  Appendix  C),  the  view  factor  for  a disk  of  radius  R,  at  a point  a distance  H along  the 
axis  normal  to  the  disk,  is 


Hence  Q = 0.592.  We  also  assume  that  T^  = 25°C  = 298  K,  and  that  = 0.85. 
Eqs(Dl)  and  (D3)  that 


(D4) 

It  follows  from 


T4  _ 


(1-0)7/ 


(D5) 


If  F = 1,  this  yields  T^j  = 493 °C,  while  F = 0.9  yields  T^  = 472.5 °C.  The  device  surface  cannot  be 
much  colder  than  500°C,  since  it  is  observed  to  glow  red.  If  T^j  = 472.5 °C,  then 

s 9 kWW. 


We  will  assume  that  the  factor  F remains  constant  for  all  irradiations.  Thus  for  the  case  = 25 
kW/m“,  <^rad*  “ kW/m",  and  Eq.(D5)  implies  that  Tj  = 840.3  K = 567. 1°C. 

We  proceed  the  same  way  for  all  four  cases,  and  obtain  the  values  in  the  first  six  columns  of  Table  Dl 
(with  fluxes  given  in  kW/m").  Column  6,  marked  <t>^,  is  the  convective  flux,  found  as  the  difference 
between  and  These  values  of  <t>^  are  plotted  ys  and  a smooth  curve  passed  through  these 
points,  including  the  point  at  the  origin.  That  yields  the  smoothed  convective  fluxes  (f)^,  given  in  Col. 7. 
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Table  Dl.  Measured  and  inferred  quantities  from  the  ignition  experiments. 


^lot 

</>rad/F 

Td(K) 

Td(°C) 

^rad 

<i>c 

<f>c 

hf 

Tg(K) 

h 

f 

18 

10 

745.7 

412.5 

9.0 

9.0 

8.9 

19.67 

698 

22.0 

0.893 

25 

16 

840.3 

567.1 

14.4 

10.6 

10.7 

19.56 

760 

22.9 

0.854 

34.5 

23 

921.0 

647.8 

20.7 

13.2 

13.0 

20.71 

837 

23.9 

0.865 

44 

32 

1000.9 

727.7 

28.8 

15.2 

15.2 

21.48 

907 

24.8 

0.866 

The  convective  flux  can  be  written  in  the  form 


(D6) 


where  Tg  is  the  purge  gas  temperature,  and  T^.  is  the  (cold)  gauge  temperature.  In  order  to  make 
progress  we  make  one  further  plausible  assumption;  the  purging  gas  takes  up  a (constant)  fraction  of  the 
total  energy  delivered  to  the  device.  We  can  express  this  as 

6,  = /e.  ff>7) 


where  0 is  the  temperature,  referred  to  the  gauge  temperature: 


Thus,  eqs  (D6)  and  (D7)  yield 


e,  = 7,-r, 


4>,  = hQ^  =fhe. 


(D8) 

(D9) 


Since  we  have  T^  (col.  4)  and  0^  (col.  7),  we  readily  infer  the  factor  fh  from  Eq.(D9).  This  is  given  as 
col.  8. 


Next,  we  must  find  the  heat  transfer  coefficient  h for  stagnant  flow.  The  purge  air  comes  down,  strikes 
the  fabric,  and  must  move  away  radially.  This  configuration  is  approximated  by  the  standard  problem 
of  stagnant  flow.  The  heat  transfer  coefficient  for  this  case  is  given  in  Kacag  et  al  (ref.  [26]) 
ref  [26]. 


Their  Eq.(2.176),  on  page  2-59,  gives  the  Nusselt  number: 


' Nu  ' 


0.767 


V-cPe 


\0.43 


(DIO) 


which  was  found  by  Cohen  (ref  [27]).  The  subscripts  w and  e,  above,  stand  for  "wall"  and  "edge"  (of 
the  boundary  layer),  respectively.  The  second  factor  on  the  right-hand  side  is  approximately  one,  in  this 
case.  The  Reynolds  number  is 
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Re  = 


(Dll) 


u i. 

•O  f 


where  is  the  characteristic  length.  One  might  think  that  the  characteristic  length  here  is  the  diameter 
of  the  opening.  In  fact,  however,  Cohen  [27]  gives  Re  in  the  form 


Re  = — ^ — 
dx  V 


(D12) 


Since  the  (vertical)  velocity  goes  from  to  zero  in  the  distance  6,  we  may  write 


du^ 

'dx  ^ 'b' 


(D13) 


then  Eq.(D12)  indicates  that  the  proper  to  use  here  is  6.  u„  is  readily  inferred  from  the  volumetric 
flow,  dV/dt  = 0.5  liters/min.  Thus 


Ti  r^uD^  lid^uT 

(D14) 

T~a  g 

which  yields 

= 6.28(7^/7^)  cm/sec 

(D15) 

The  Nusselt  number  is  given  by 

SI  hb 

Nu  = — 

(D16) 

K 

where  6 is  the  stand-off  distance  between  the  heater  and  the  substrate,  5.4  mm.  Combining  these,  and 
knowing  jt(T)  and  ^(T)  for  air,  we  obtain  h(T  ),  as  shown  in  Fig.  17.  Note  that  this  is  independent  of  any 
estimates  of  F,  the  radiative  and  convective  fluxes,  etc. 

We  now  proceed  as  follows:  we  guess  a value  for  Tg;  corresponding  to  this  we  have  h(Tg),  and  we  then 

find  h(Tg)(Tg  - TJ  s <^*.  We  must  do  this  until  </>*  = <l)^,  as  given  in  Table  Dl.  This  procedure 
converges  quite  rapidly,  and  we  find  the  values  of  Tg  and  h(Tg)  shown  in  columns  9 and  10.  Finally, 
dividing  hf  by  the  now-known  values  of  h,  we  obtain  the  values  of  f shown  in  col.  11.  There  are  two 
things  to  be  noted  about  f:  it  is  almost  as  large  as  it  can  get  (i.e.,  close  to  unity),  and  it  is  approximately 
independent  of 
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APPENDIX  E 

Thermophysical  Data  for  Cotton 

The  fabric  is  a weave  of  cotton  fibers.  Cotton  itself  is  a flexible,  hollow  tube  of  cellulose;  the  central 
channel  is  called  the  lumen,  and  occupies  a fraction  t]  of  the  total  volume.  Thus  we  would  expect  that 
the  density  of  cotton  is  about 

p(cotton)  = (1-t;)p3, 

where  is  the  density  of  the  solid  (largely  a-cellulose).  Measurements  yield  rj  = 0.2  - 0.4.  However, 
refllT]  gives  the  following  data  on  page  V-122  (the  references  given  there  to  the  original  authors  are 
omitted  here  for  brevity): 


Material  Density  (g/mf) 


Cellulose  I 
Cellulose  II 
Cellulose  III 
Cotton 


1.582- 1.630 

1.583- 1.62 
1.61 

1.545-1.585 


Thus,  the  density  of  cotton  appears  to  be  very  nearly  the  same  as  that  of  the  solid  (cellulose)  implying 
that  t;  » 0.  We  will  hereafter  ignore  the  apparently  small  difference  between  a-cellulose  and  cotton,  and 
take  the  (mean)  density  to  be  p^  = 1.565  ± 0.02  g/mf. 

For  many  polymers,  there  is  only  a weak  dependence  of  the  thermal  conductivity  on  T,  between  100  and 
300  K (see  Fig. 68b,  ref  [25]).  One  can  get  an  idea  of  the  variability  of  the  thermal  conductivity  of  cotton 
and  of  cotton  fabrics  from  Figs.AA-AC,  in  ref  [26]. 

Measurements  at  NIST  by  J.R.  Lawson  (private  communication)  have  shown  that  the  density  of  #10  duck 
is  p = 0.6  g/cm^,  that  of  #6  duck  is  0.72  g/cm^,  and  that  of  all  the  other  cotton  duck  fabrics  measured 
is 

Pf  = p(fabric)  = 0.62  g/cm^  = 620  kg/m^. 

The  "void  fraction"  of  the  fabric  is  With  Pj  = 1.565,  Pf  = 0.62,  and  with  p^  = density  of  air  at 
STP  = 1.774x10'^  g/cm^,  the  relationship 

Pf  = (1  - <J>)p3  -I-  ^p^ 


yields  $ = 0.6045. 

Next,  consider  the  specific  heat.  Again,  the  values  given  in  ref  [17]  are: 
Material  Specific  Heat  (J/g-K) 


Cellulose 

Cotton 


1.34 

1.22 

1.15-1.18  (0-34°C) 
1.32  (0-100°C) 
1.327-1.251 
1.273-1.35 
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Where  the  temperature  range  is  not  listed,  it  is  assumed  to  be  the  ambient  (20  or  25°C).  We  will  thus 
not  be  far  wrong  it  we  take 

C3=1.3J/g-K, 

where  the  subscript  "s"  stands  for  "solid"  (cotton  or  a-cellulose).  We  expect  that,  just  as  was  the  case 
for  density, 

Cf  = c(fabric)  = #0^;^  + (1  - ^jc^. 


Hence 

Cf  = 1.122  J/g-K,  atT«  300K. 

This  author  has  not  been  able  to  find  the  temperature  dependence  Cs(T);  we  will  defer  that  question  for 
the  moment. 

Finally,  we  come  to  the  thermal  conductivity,  k.  This  is  a quantity  which  is  notoriously  difficult  to  obtain 
accurately.  Figure  18,  from  ref  [26],  shows  the  large  variations  in  thermal  conductivity  depending  on 
measuring  conditions.  Many  of  the  low  values  (nos.  5,  6,  and  8,  for  example)  were  measured  in  a 
vacuum.  Numbers  4 and  9,  on  the  other  hand,  were  similar  specimens,  measured  in  air  at  25 °C; 
however,  #9  had  over  three  times  the  density  of  #4  (and  would  therefore  have  been  expected  to  have  a 
substantially  higher  thermal  conductivity!). 

We  list  here  four  sources  for  the  thermal  conductivity: 

a.  The  apparently  most  consistent  data  from  [26]  --  curves  1,  2,  and  3,  and  point  7,  give 

x = 0.0365  W/m-K  at  T = 306  K = 33‘’C 

b.  We  have  the  following  data  for  cotton  paper,  from  ref.  [17]: 

T (°C)  X (W/m-K) 


30  0.049 

40  0.071 

50  0.084 

60  0.090 

A curve-fit  to  these  points,  if  extrapolated,  would  predict  that  x vanishes  at  (and  below)  T = 10°C, 
which  is  nonsense. 

c.  Ref  [27]  gives: 

T (°C)  X (W/m-K) 

— The  density  of  the  material  for  these 

0 0.056  measurements  is  listed  as  81  kg/m^;  hence 

20  0.058  $ = 0.9493.  Then  Eq.(50)  gives 

100  0.067 

X = 0.03408X,  -b  0.9828X  P) 
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while  the  tabular  data  are  reasonably  well  fitted  by 
/c(T)  = 0.056  + l.lxia^T,  with  T in  °C 


(f2) 

We  have  k^^(25°C)  = 0.02624  W/m-K;  now  if  the  data  in  the  table  are  the  conductivities  for  the  light 
cotton/air  mixture,  then  at  20°C,  with  Eq.(El)  yields  = 0.945,  which  is  an  unreasonably 

large  value.  According  to  Kunii  (ref  [4]),  the  value  for  the  gas-phase  thermal  conductivity  to  be  used  in 
Eq.(50)  is  1.7  times  what  it  is  in  the  ambient.  If  we  use  = 1.7x0.02624  = 0.04461,  the  resulting 
value  of  is  = 0.415  W/m-K,  still  very  large.  If,  on  the  other  hand,  the  listed  values  are  for  k^,  then 
Xg  = 0.058  and  k„  = imply  that  k = 0.0278,  while  for  = 1.7>c^,  k = 0.0458.  Thus  none  of  the 
combinations  is  plausible. 

d.  T.J.  Ohlemiller  has  made  measurements  of  the  thermal  conductivity  of  the  (#12)  cotton  duck;  the 
apparatus  only  works  properly  when  the  sample  is  thermally  thick,  however.  Therefore  he  carried  out 
a series  of  measurements,  with  an  increasing  number  of  layers  of  the  fabric.  It  was  then  possible  to  infer 
the  asymptotic  value  which  would  be  reached  if  the  number  of  layers  had  been  increased  to  a very  large 
number:  It  was  assumed  that 

x(n)  = x„[l  - exp(-ne)], 

where  0 = dimensionless  thickness  of  one  layer,  and  n = number  of  layers.  The  data  were 


n K 


0 0 
6 0.346 

12  0.533 

18  0.683 


The  thermal  conductivity  is  given  here 
in  units  of  BTU-in/ft”-hr-°F.  The 
asymptotic  value  of  the  series  at  the 
left  is  0.9.  Converted  to  units  of 
J/m-K,  that  is  = 0.13  J/m-K. 


X = 0.13  is  about  twice  the  value  (0.056)  found  in  the  Handbook,  as  quoted  in  c.,  just  above.  Recall 
that  for  this  fabric,  the  void  fraction  is  0.6045.  With  $ = 0.6045,  Eq.(50)  gives 


X = 0.28505/Cg  -h  0.84554xg^ 

If  we  use  /Cg^  = x^jr,  we  would  then  infer  that  Xg  = 0.378,  which  is  an  order  of  magnitude  greater  than 
the  earlier  estimates.  If  we  use  = 1.7xair,  then  x = 0.13  implies  Xg  = 0.324,  only  a little  smaller. 
The  latter  value  is  not  too  different  from  the  0.415  found  in  section  c.,  above. 


We  thus  take  x^,  = XgCT^  = 25 °C)  « 0.324  W/m-K 

and  x(TJ  = 0.13  W/m-K 

Assuming  that  the  temperature  dependence  is  like  that  given  by  Eq.(E2),  we  find  that 

x(T)  = 0.125  -1-  2.457x10-^  T,  (E3) 


with  T in  "C.  Alternatively,  we  might  use  Eq.(50)  in  its  more  general  form  (that  is,  include  the 
temperature  dependence): 

x(T)  = 0.28505xg(T)  -I-  0.84554(1.7)Xairn')  (E4) 
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For  air,  it(T)  is  given  by  Eq(52)  (see  ref  [28]).  For  the  solid,  a reasonable  assumption  is  that  the 
conductivity  is  proportional  to  the  absolute  temperature,  and  that  it  is  also  proportional  to  the  density. 
Thus  Ks(T)  is  given  by 

K,(J)  = pT/pJ,  = 0.3237  pT/pJ,  (E5) 

We  now  turn  to  the  question  of  the  temperature  dependence  of  Cf.  For  a number  of  materials,  the 
thermal  diffusivity  is  relatively  insensitive  to  T.  This  is,  in  particular,  the  case  for  wood  (see  [15],  [16]) 
which  consists,  to  a significant  extent,  of  cellulose  and  related  compounds.  If  we  assume  that  to  be  the 
case  for  cotton,  then  we  could  write 

ic(T)  = apcCT);  (E6) 

c(T)  is  generally  much  easier  to  measure  than  /<(T),  and  this  would  be  a relatively  good  way  to  obtain 
ic(T);  the  irony  is  that  we  do  not  have  c(T)  for  cotton  or  cellulose. 

If  we  take  Kf  = 0.13  for  the  fabric  at  ambient  temperature,  we  now  find  that 

a = K/pc  * 0.13/(620)(1122)  = 1.87x10“^  m^/s 

We  therefore  have,  finally, 

c(T)  = K/pa  * x(T)/620(  1.87x10-'^)  = 8630/t(T)  J/kg-K 
where  x(T)  is  given  by  Eq(E3)  (or  by  (E4),  if  preferred). 
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